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1 .  INTRODUCTION 


.1.  Methods  for  Reduced  Order  Modeling  with  Applications  to  Power  Systems 


An  object  of  study  in  this  thesis  is  the  relation  between 
singular  perturbation,  modal  methods  and  coherency  when  used  for  modal 
order  reduction  of  electromechanical  models  of  power  systems.  We  will 
show  that  all  these  methods  are  closely  related  if  coherency  is  suitably 


defined . 


Modal  method,  originally  proposed  by  Davison  [1]  and  further 


analyzed  by  many  researchers  [2,3]  consists  of  approximating  the  system 

x  =  Ax  ,  dim  A  =  n  (1-1 

by  one  of  lower  order 

z  =  Fz  ,  dim  F  =  r .  (1.2 

It  can  be  shown  [4,5]  that  modal  methods  can  be  understood  as  a  particular 
case  of  aggregation  [6,7],  i.e.  there  exists  a  matrix  K,  such  that 

z  =  Kx  (1-3 


FK  =  KA. 


Aggregation  matrix  F  is  then  given  as 


F  =  KAK+  ,  KK+  *  I 


In  [5]  it  is  shown  that  K  has  to  be  of  the  form 


K  =  M(I  0)X 


where  M  is  an  arbitrary  nonsingular  rxr  matrix  and  X  is  the  nxn 
matrix  of  right  eigenvectors  of  A.  With  such  K,  and  <?(•)  denoting 


2 


spectrum  of  a  matrix. 


<J(F)  <=  cr  (A)  . 


In  modal  methods  of  [1]  M  =  Xrr,  which  is  r  -  dimensions  1  principal 


(1-7) 


minor  of  X  ordered  so  that  the  first  r  columns  correspond  to  the  eigen¬ 


values  to  be  retained  in  F,  and 


F  A11  +  A12X21X11 


(1.8) 


where  X,,  =  X  and  X„ ,  is  "21"  block  in  the  matrix  X. 

11  rr  21 

In  the  singular  perturbation  approach  [8],  it  is  assumed  that 
the  state  vector  x  consists  of  subvector  x^,  with  r  basically  slow 
variables,  and  of  the  complementary  vector  x^  having  variables  with 
mixed  fast  and  slow  dynamics.  The  presense  of  different  speeds  of 
response  in  these  two  sets  of  variables  is  indicated  by  writing  the  model 
(1.1)  in  the  form 


*1  '  A11X1  +  A12X2 


S  "  *21*1  +  A22X2 

where  £  >  0  is  a  given  small  number.  The  slow  behavior  of  the  system 
is  analyzed  via  the  descriptor  variable  system  [9,10] 


(1.9) 


x  =  A,,x  +  A  xc 
s  11  s  12  f 


0  A21xs  +  A22xf 


(1.10) 


obtained  from  (1.9)  by  setting  £  to  zero.  Assuming  A2?  exists  the  slow 


dynamics  is  modeled  as 


xs  ^All  "  A12  A22  A2  l^xs 


(1-11) 


The  slow  eigenvalues  and  the  corresponding  slow  elgenspace  of  (1.9)  are 
continuous  in  £  [10].  Consequently  the  slow  component  of  system  responses 
is  also  continuous  in  £  [8].  Therefore  the  slow  behavior  of  (1.9)  is 
used  as  an  approximation  of  the  slow  behavior  of  (1.9),  whenever  £  is 
sufficiently  small.  If 


then 

X1  =  Xs  + 

X2  =  Xs  +  Xf  ‘  Xs(0)  +  °(£)  (1> 1 

In  the  model  (1.9)  there  exists  not  only  an  eigenvalue  separation 
(going  to  infinity  as  £  tends  to  zero)  [8],  but  also  there  is  a 
difference  in  speed  of  variables  in  x^  and  x^ •  Such  systems  are  called 
[11]  state  separable  systems.  Systems  which  have  eigenvalue  separation 
between  slow  and  fast  eigenvalues,  but  state  vector  can  not  be  separated 
into  x^  and  based  on  different  speeds  of  response  are  called  [11] 
mixed  state  systems.  For  such  systems  and  systems  with  larger  £  an 
iterative  separation  of  time  scales  of  [12,13]  can  be  used. 

The  problem  of  the  time  scale  separation  consists  of  finding 
matrices  L  and  H  in  the  state  transformation  matrix 


(1.1 


such  that  L  satisfies  a  generalized  matrix  Riccati  equation 


4 


is 


R(L)  =  A22L  -  IAn-  LA12L+  A21  =  0 


and  H  satisfies  a  Lyapunov  type  equation 


P  (H)  =  -BjH  +  HB2  +  A12  -  0 


In  the  above  equations  and  B2  are  defined  via 


(1.15) 


(1.16) 


TAT 


-1 


B1  ° 


0  B0 


(1.17) 


The  solutions  L  and  H  have  to  be  such  that 


Imb^I  <  |  X(B2)|  . 


(1.18) 


A  central  problem  in  so  defined  time  scale  decomposition  is  existence  and 
form  of  the  solution  for  L.  This  problem  has  been  studied  in  [31,14,15,13] 
It  has  been  shown  that  the  generalized,  like  the  standard  Riccati  equation, 
has  many  solutions;  furthermore,  whenever  there  is  a  solution  for  L  the 
existence  of  a  solution  for  H  is  guaranteed.  Several  numerical  schemes 
have  been  proposed  for  solving  the  Riccati  equation  [12,13].  The 
typical  sequence  has  the  form 


-1 


Lk+1  (A22+LA12)  R^)  +  1^ 


(1.19) 


Both  modal  methods  [16,17,18]  and  the  singular  perturbation 
methods  [11]  have  been  applied  to  power  system  problems.  Independently 
of  these  general  system  concepts,  in  the  analysis  of  power  systems  a 
specific  method  called  coherency  has  been  used  for  model  order  reduction. 


5 


This  concept  relies  on  an  empirically  observed  phenomenon  that  after  a 
disturbance  system  responds  in  groups  of  machines,  such  that  in  each  group 
all  the  machines  have  the  same  response.  These  groups  are  called  coherent 
groups.  A  numerical  algorithm  for  identifying  coherent  groups  based  on 
simulation  of  machine  responses  after  a  disturbance  and  checking  which 
machines  satisfy  the  above  property  is  given  in  [19].  In  [20,21]  it  is 
shown  that  coherency  of  a  group  of  machines  is  equivalent  to  the  reachable 
subspace  of  the  disturbance  being  in  the  null  space  of  the  output  matrix 
defined  to  give  as  output  the  differences  of  machine  angles  in  the  groups 
In  [16]  it  is  indicated  that  coherency  can  be  identified  in  a  special 
case  of  the  modal  method  analysis,  when  of  (1.3)  has  a  special 

structure,  namely  if  each  row  has  all  elements  zero  except  for  one  which 
is  1.  Other  approaches  to  coherency  include  [22,23];  an  extensive  list  of 
references  is  contained  in  [24] . 

Having  established  the  connection  between  modal  methods,  singular 
perturbation  and  iterative  time  scale  decomposition  in  [14],  it  is  the 
above  result  of  [16]  that  motivated  the  work  on  defining  partial 
coherency  in  the  framework  of  the  time  scale  decomposition  [25],  i.e.  in 
terms  of  the  Riccati  equation.  The  alternative  formulation  of  the 
coherency  led  to  a  number  of  new  results  surveyed  in  the  next  section. 

1.2.  Contributions  of  the  Thesis 

The  major  contribution  of  this  thesis  is  in  establishing  a  link 
between  time  scales,  coherency  and  coupling  between  areas.  It  is  shown 
that  if  coherency  is  defined  with  respect  to  the  slow  modes  of  the  system, 
resulting  areas  are  weakly  coupled.  In  this  case  interarea  motion 


represents  the  slow  dynamics  of  the  system  and  the  intermachine  oscillations 
within  coherent  groups  represent  systems  fast  dynamics.  For  the  first 
time  it  is  argued  analytically  and  then  confirmed  experimentally  that  so 
defined  areas  are  relatively  robust  with  respect  to  parameter  changes  and 
fault  location  in  the  system.  The  consideration  of  coherency  is  completed 
by  presenting  a  novel  algorithm  for  identifying  coherent  areas  which 
beneficially  exploits  an  analytical  definition  of  coherency  via  the  Riccati 
equation. 

This  work  has  been  motivated  by  desire  to  overcome  limitations 
of  existing  methods  for  reduced  order  modeling  of  power  systems.  For 
example,  all  existing  methods  for  identifying  coherency  [20,21,19]  depend 
heavily  on  the  fault  location,  making  the  identified  areas  valid  only 
for  the  particular  disturbance.  The  analytical  work  of  [20,21]  qualifies 
coherency  analytically,  but  suffers  from  the  same  fault  dependance  and 
lacks  an  efficient  algorithm  for  identifying  areas.  The  modal  method  of 
[18]  is  based  on  the  linearized  model  of  a  part  of  the  system  (the  one 
outside  of  a  specified  subsystem,  usually  called  the  study  area).  As  such 
the  approach  can  give  errorneous  results  in  the  case  of  stronger  coupling 
between  the  two  parts.  The  other  common  objection  of  the  method  is  that 
the  reduced  order  system  looses  any  physical  meaning.  In  the  modal 
approach  of  [16]  the  first  objection  is  eliminated  by  considering 
linearized  model  of  the  whole  system.  Furthermore,  the  connection  with 
coherency  is  indicated  in  case  of  the  matrix  X21X11^  hav*n§  the  special 
structure,  as  discussed  earlier.  However,  neither  explicit  characteriza¬ 
tion  of  this  coherency  nor  the  method  of  manipulating  the  system  so  as  to 
achieve  this  desired  form  of  L  given.  In  [26]  coherency  is  defined  as 


partial  coherency  with  respect  to  a  desired  subset  of  eigenvalues. 

The  algorithm  for  identifying  coherent  areas  and  the  properties  of  those 
are  not  elaborated  on.  Furthermore,  none  of  the  known  works  analyses 
analytically  robustness  of  the  area  boundaries  with  respect  to  parameter 
and  load  changes . 


1.3.  Chapter  Preview 

Chapter  2  establishes  a  background  for  the  analysis  of  coherency 
in  power  systems,  which  is  the  main  theme  of  this  study.  Ihis  background 
is  found  in  the  spectrum  separation  problem  formulated  in  terms  of  a  - 
general  matrix  Riccati  equation.  Properties  of  its  solution  are  studied 
and  as  a  result  a  relation  between  singular  perturbations,  iterative 
separation  of  time  scales  and  the  modal  methods  is  established.  An 
algorithm  for  the  time  scale  decomposition  of  weakly  nonlinear  systems  is 
proposed . 

A  short  overview  of  the  existing  methods  for  the  coherency 
analysis  of  power  systems  is  given  in  Chapter  3.  Then  a  new  defini¬ 
tion  of  coherency  as  a  partial  coherency  with  respe-'t  to  a  given 
subspectrum  is  given  in  the  framework  of  the  spectrum  separation  problem. 

J 

A  class  of  systems  called  r-decomposable  systems  is  introduced  to  study 
the  idealized  coherency.  The  gained  insight  is  then  used  in  constructing 
an  algorithm  for  coherency  identification  of  realistic  (large)  power 
systems.  An  algebraic  criterion  is  given  for  checking  the  validity  of  the 
use  of  so  called  "area"  variables,  frequently  used  by  power  system  analysts. 


Aggregation  of  linear  and  nonlinear  power  system  models  based 
on  the  partial  coherency  is  studied  in  Chapter  4.  When,  instead  of  an 
arbitrary  partial  coherency,  the  so  called  slow  coherency  is  used,  the 
resulting  areas  are  insensitive  to  small  variations  in  network  parameters. 
Furthermore,  in  this  case  the  areas  are  weakly  coupled.  Writing  the 
electromechanical  model  in  terms  of  the  area  variables  and  machine  angle 
differences  (the  differences  of  angles  within  coherent  areas)  transforms 
the  originally  mixed  state  electromechanical  model  into  the  singularly 
perturbed  form.  Use  of  the  singular  perturbation  theory  and  the  weak 
coupling  for  the  simulation  of  the  systems  responses  is  illustrated  on  a 
nontrivial  power  system  example. 

Further  uses  of  coherency  in  the  direct  transient  stability 
analysis  and  other  possibilities  of  extending  this  study  are  discussed 
in  Chapter  5. 
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2.  TIME  SCALE  DECOMPOSITION  VIA  RICCATI  APPROACH 
2.1.  Problem  Formulation 

In  this  chapter  a  group  of  results  will  be  given  concerning 
linear  systems,  among  which  we  will  have:  the  form  and  the  properties  of 
solutions  of  Riccati  equation  (1.15);  a  relation  between  the  time  scale 
decomposition,  modal  methods  and  singular  perturbation  methods  for  reduced 
order  modeling;  and  convergence  properties  of  a  class  of  iterative 
algorithms  (1.19).  For  weakly  nonlinear  systems  an  algorithm  for  simulation 
is  proposed,  which  takes  advantage  of  time  scale  decomposition  of  the 
linear  part.  The  algorithm  is  related  on  an  example  system  with  an 
alternative,  purely  numerical  algorithm. 


c 


2.2.  Time  Scales  and  Generalized  Riccati  Equation 

We  will  first  define  a  general  spectrum  separation  problem 

and  then  as  a  special  case  treat  the  time  scale  decomposition.  We  start 

with  the  spectrum  separation  problem  as  given  in  [27]  and  then  proceed 

with  physically  more  meaningful  formulation  of  the  same  problem  via  solution 

of  a  Riccati  equation.  The  problem  of  spectrum  separation  can  be  defined 

as  one  of  finding  two  lower  order  matrices  and  F^  such  that  F^ 

reproduces  a  given  subspectrum  of  A,  cr  ,  and  F^  reproduces  the  complementary 
c  c 

subspectrum  o^,  ar  U  a  =<  o'(A).  We  further  want: 

Assumption  2.1:  H  =  0. 

According  to  [27,  first  decomposition  theorem],  the  problem  of 

spectrum  separation  is  always  solvable  under  the  Assumption  2.1.  To  each 

of  the  subspectra  there  corresponds  a  unique  A-invariant  subspace, Sf  c  X 

c  . 

of  dimension  r  and  S  -  X  of  dimension  n  -  r,  where  X  is  n-dimensiona 1 

r 


\ 
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( 

p 

I 


i 


vector  space  ,  such  that  the  spectrum  of  A  restricted  to  S^.  is  a  ,  and  the 

c  c 

spectrum  of  A  restricted  to  S is  c  .  Then  A  in  this  new  basis  is  of  the 


form: 


A  = 
r 


0 


(2.1) 


where  c(F^)  =  cr^,  ct(F2)  =  O’  •  To  verify  above  form  let  us  note  that  for 

all  column  vectors  in  V,  we  have  AV,€  S  •  Hence,  there  exists  a  matrix 

1  1  r 

such  that 

AV1  =  V1F1  »  (2>2) 

but  this  defines  eigenvalue-eigenvector  problem  up  to  a  similarity 
transformation  =  X^K^,  where  X^  is  the  matrix  of  r  eigenvectors  corre- 

Q 

sponding  to  Similarly,  for  all  vectors  in  ,  AV^  and  there 

exists  a  matrix  F2  such  that 

AV2  =  V2F2  (2.3) 

Q 

with  a(F2)  =  cr^  and  V2  =  ^2^  for  eigenvectors  X7  and  some  nonsingular  . 

The  relation  between  the  states  of  equivalent  systems  in  which 
A  and  Av  are  state  matrices  is  given  by 

x  =  Vxv  .  (2.4) 

Notice  that  in  order  to  achieve  the  spectrum  separation  it  is  sufficient 
to  consider  only  the  invariant  subspace  corresponding  to  the  given 
subspectrum  <j  ,  and  complement  the  space  by  any  independent  n-r-dimensional 
subspace.  Of  course,  in  this  case  one  gets  only  a  blocktriangular  form  of  A  . 

■A1 

S  c  X  is  A-inva riant  if  AS  c  S 
r  r  r 
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Let  us  consider  now  the  transformation  of  (1.1)  in  the  form 


x  =  [Vx  V2]xy 


(2.5) 


where  ImV^  is  A-invariant  r-dimensional  subspace  of  Z  corresponding 
to  and  Im  V2  spans  n-r  dimensional  subspace  of  Z  such  that  it  does 
not  contain  any  of  the  vectors  in  the  subspace  ImV^,  but  otherwise 
arbitrary. 


Remark  2.2:  Since  has  full  rank,  there  always  exists  an  ordering  of 

M 

states  in  x  for  which  V.,  of  V..  =  is  nonsingular.  Due  to  this 

ii  1  V 

_ 21 

remark  we  assume  that  x  in  (2.5)  is  ordered  so  that  is  nonsingular. 


That  allows  a  particularly  simple  ^2*  ^2  =  [ij  ‘  choice  of 

and  V2  the  transformation  (2.5)  can  be  rewritten  as 


V  V_1  I 

21  11  j 


°  Kl  ° 


(2.6) 


The  transformed  state  matrix  is 


0  I  -L  I 


0  I  I  0  I  0  V--  0  F,  A,» 

A  11  =  (2.7) 


L  I 


J  L°  L°  F2  . 


where 


L  S  V21Vli1 


(2-8) 


F1  ~  V11  (A11 +A12L)V11  V11  B1V11 


(2.9) 


F 2  ”  L  A 12 


(2 . 10) 


12 


It  follows  from  (2.7)  that  the  spectrum  separation  can  be  achieved  with 
the  transformation 


x  = 


I  0 
L  I 


x  =  T.  ^x 
v  L  v 


(2.11) 


with  the  only  requirement  that  Im 


0 


is  A-invariant  corresponding  to 


a  .  The  significance  of  this  transformation  is  that  the  states  of  x 
r  v 

retain  the  physical  meaning  if  the  spectrum  separation  problem  results 
in  L  =  0,  i.e.  if  the  original  system  consists  of  two  decoupled 
(through  A^)  subsystems.  In  cases  of  nonzero  coupling  between  the 
subsystems,  the  first  r  states  of  xy  still  retain  their  meaning,  while 
only  the  last  n-r  states  are  changed.  After  this  transformation  is 
applied,  we  get 


A  = 
v 


B, 


12 


(2.12) 


where  and  B2  are  defined  by  (2.9)  and  (2.10). 

From  (2.7)  it  also  follows  that  in  solving  the  spectrum  separa¬ 
tion  problem  it  is  necessary  to  make  "21"  block  in  the  transformed  matrix 
equal  zero.  By  writing  the  equation  for  "21"  block  of  (2.7)  we  get 


R(L)  ’  -  LAn  +  A22L  -  LA12L  +  A21  =  0  (2.13) 

This  equation  has  the  form  of  a  generalized  Riccati  equation.  In  the 

T 

case  when  =  -  A22  it  becomes  the  standard  algebraic  Riccati 
equation  found  in  optimal  control  and  estimation  theory.  For  later 
reference  we  write  this  equation  in  the  form  [28,29] 
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0 


R  (L)  =  -  LF  -  FTL  +  LBR-1BTL  -  Q 
s 

for  which  underlining  problem  is  one  of  finding  control  u  of  the 
system 

x  =  Ax  +  Bu 

so  as  to  minimize  J(x,u)  given  by 

CO 

J(x,u)  =  J*  (xTQx  +  uTRu)dt 
o 

The  matrix  E  associated  with  (2.14),  the  same  way  A  is  associated  to 
(2.13),  has  the  form 


(2.14) 


(2.15) 


(2.16) 


E  = 


-1  T 

A  -BR  B 


-Q  -A 


(2.17) 


Therefore  equation  (2.14)  is  a  special  case  in  the  spectrum  separation 
problem.  Several  standard  results  of  (2.14)  will  be  rephrased  in  the 
light  of  the  results  on  existence, form  and  uniqueness  of  the  solution 
of  generalized  Riccati  equation  (2.13),  which  are  coming  next. 

Theorem  2.3: 

Let  an  r-element  subspectrum  of  A,  be  given,  with  no  common 

c  c  c 

eigenvalues  with  its  complement  a  ,  fl  7s  0,  U  =  cr(A)  .  Let 

be  the  corresponding  A-invariant  subspace  of  dimension  r  with  its  basis 

nx  *■ 

given  as  columns  of  a  matrix  V€R  ",  i.e. 


AV  =  VF, 


(2.18) 


implies 


a<Fl>  *  c c 


(2.19) 


(i)  the  solution  L  of  the  Riccati  equation  (2.13),  such  that 

rv 

ct(bi)  =  cr  exist  if  and  only  if  the  matrix  V.  of  V  =  is  nonsingular. 

r  Lv2  J 

This  solution  is  then  given  as 


(2.20 


(ii)  When  columns  of  V  are  eigenvectors  of  A,  then  columns  of  V 


are  eigenvectors  of  B^. 


Proof:  Equation  (1.13)  can  be  rewritten  as 


R(L)  =  [-L  I] A 


El- 


(2.21 


The  equation  will  be  satisfied  if  and  only  if 


(a)  Im  C  Ker  A  ,  or 
L 


(b)  Im  A  ^  c  Ker  [-L  I]  =  Im 

Li 


But  since  (a)  is  already  included  in  (b)  we  proceed  with  (b) 


’ll  .  T  V 
[ L J  ~  A  In  [l_ 


c  Ker  !-L  I]  =  Im 


UJ  * 


(2.22; 


i.e.  for  any  solution  L,  Im  must  be  A-invariant.  We  now  want  to 

LL  J 

relate  this  subspace  to  the  spectrum  a  .  More  specifically  we  will 
be  given  and  the  corresponding  A-invariant  subspace  S  with  a  basis 


V£r“*l  ,  and  we  want  to  see  if  there  exists  an  L  for  which  1"  will  span 

Li 

Sr .  Clearly,  if  and  only  if  is  nonsingular  such  an  L  exists,  since 


from  the  requirement  that  V  and  t  span  the  same  subspace,  there  exists 


a  matrix  K€R  for  which 


(2.23 


From  rank  1=  r=  rank  V^K  it  follows  that  both  and  K  have  to  be  of  full 


.-1 


rank.  Then  K = ;  which  proves  (2.20).  To  prove  a(B^) =  we  write  as 


Ti 


Bl=  [I  0]  A  j^Lj 


(2.2  4 


Due  to  (2.18)  and  (2.23)  with  K  =  ,  we  have 


B1  =  [I  0]  AVV^1  =  [I  0]  VF1V^1  =  V][F1V^1 


which  completes  the  proof  of  (i)  and  (ii). 

The  proof  given  here  is  more  general  than  those  given  in 
[31,15,14].  In  [31]  only  the  case  of  diagonalizable  matrix  A  has  been 
treated  and  a  particular , eigenvector  basis  for  was  assumed.  In 
[14,15]  sufficiency  has  been  proved.  Namely  if  V  spans  an  invariant 
subspace  then  the  solution  L  is  of  the  form  (2.20).  Here  it  is  proved 
that  for  L  to  exist  it  is  also  necessary  that  V  spans  an  invariant 
subspace  of  A.  The  form  of  the  solution  also  applies  to  the  standard 
Riccati  equation  (2.14)  with  E  of  (2.17)  playing  the  role  of  A. 

A  different  proof  for  this  equation  is  given  in  [30] . 

Remark  2.4:  From  the  geometrical  arguments  used  in  the  proof  of  the 
theorem,  it  should  be  clear  that  L  does  not  depend  on  the  particular 
choice  of  basis  for  an  invariant  subspace  corresponding  to  a  given  spectrum 
We  will  however  demonstrate  this  fact  in  matrix  terms,  which  are  more 
suitable  for  computational  purposes. 

Proposition  2.5:  Each  matrix  L  solving  the  Riccati  equation  (2.7)  is 
independent  on  the  choice  of  V  spaning  A-invariant  subspace  corresponding 


to  a  given  spectrum. 


Proof :  Suppose  that  V  spans  A-invariant  subspace  corresponding  to  a  given 

e  ar 

Fvi' 


spectrum  cr^.  By  assumption  on  existence  of  L,  is  invertible  and 

l  =  v^;1. 


If  some  V'  ^  V  spans  the  same  subspace,  then  V'  = 


some  nonsingular  K.  In  this  case 


Lv2j 


K,  for 


L'  =  (V2K)(V1K)"1  =  L,  (2.25) 

which  completes  the  proof.  From  the  theorem  it  follows  that  the  Riccati 

equation  has  many  solutions.  In  fact  for  any  cr^  for  which  eigenvector 

basis  V  has  invertible  (due  to  the  Proposition  2.1  any  other  basis  can 

be  assumed)  there  corresponds  one  L.  For  the  time  scale  decomposition 

of  particular  interest  is  the  case  when  ct  =  o  ,  where  cr  contains  r 

r  s  s 

smallest  in  magnitude  eigenvalues  of  A.  In  this  case  the  eigenvalues 

of  A  are  divided  into  the  slow  eigenvalues  contained  in  B^  and  the  fast 

eigenvalues  contained  in  B^ •  Due  to  the  crucial  role  of  this  solution  in 

the  most  of  the  results  that  follow,  we  introduce 

Definition  2.6:  The  solution  of  Riccati  equation  (2.13)  for  which 

!  X. ( B .^ )  I  <  |X(B2)|  is  called  dichotomic  and  is  denoted  by  L^. 

Rema rk  2.7;  In  case  of  the  standard  Riccati  equation  (2.14)  and  (2.17), 

it  is  known  [29]  that  under  usual  controllability  -  observability 

assumptions,  matrix  E  has  exactly  n  eigenvalues  with  strictly  negative 

real  parts,  whose  spectrum  we  will  denote  by  a  and,  symmetric  to  those 

n 

with  respect  to  the  imaginary  axes,  n  eigenvalues  with  positive  real 
part,  whose  spectrum  will  be  Among  many  solution  of  (2.14),  of 

-IT 

particular  interest  is  the  one  for  which  cr  =  c  and  B,  =  A„  =  F  -  BR  B  L. 

r  n  1  F 

In  this  case  the  Assumption  2 . 1  is  satisfied.  By  the  Theorem  2.3  ?(A  =  cr  , 

F’  n 
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i.e.  the  closed  loop  system  matrix  is  stable.  It  is  also  easy  to  see 

that  there  exists  one  solution  L  for  which  ct(Af)  =  a*.  The  symmetry  of  L 

T  T 

can  be  concluded  from  R  (L  )  =  [R_,(L)]  =  0  and  uniqueness  of  the  solution 

5  S 

corresponding  to  the  same  spectrum. 

Remark  2.8:  Let  X  be  the  matrix  of  right  eigenvectors  (and  generalized 

I  -1 

eigenvectors)  of  A.  If  V  of  Theorem  2.1  is  V  =  ;  |  ,  then  L  =  X9 . X 1 . 

i  An  «  f  ^  *  11 

_  j  L  21 j 

and  +  ^2^21^11'  But  this  has  exactly  the  same  form  as  the 

aggregated  state  matrix  in  modal  method  approaches.  The  spectrum 
separation  problem  as  defined  by  (2.1)  and  (2.2)  indicates  one  less  known 
fact,  that  with  almost  the  same  computational  expense  needed  to  obtain 
in  modal  method  approaches  one  can  get  also  the  matrix  B^  (2.10)  which 
reproduces  the  complementary  spectrum  of  A. 

From  the  statement  of  the  theorem  and  when  Assumption  2.1  applies 
it  can  be  seen  that  the  only  condition  under  which  the  solution  of  the 
Riccati  equation  does  not  exist  is  when  is  singular.  Since  at  the  time 
when  a  spectrum  <j  is  specified  and  some  basis  of  the  corresponding 
A-invariant  subspace  is  computed  it  is  not  known  whether  will  be 

singular  or  not,  the  question  arises:  can  the  effort  spent  on  computing  V 
(which  is  considerable  as  will  be  discussed  in  chapter  3  and  4)  be  still 
used  to  achieve  the  desired  spectrum  separation,  even  if  is  singular. 

The  answer  is  yes,  due  to  the  following  Proposition  [31]. 

Proposition  2.3:  If  V  spans  r-dimensional  A-invariant  subspace  corre¬ 
sponding  to  a  given  C7^_,  then  there  always  exists  a  state  permutation 

A  A  T  A  A  /> 

x  =  Px,  A  =  PAP*  such  that  V  =  PV  is  A-invariant  and  is  invertible. 
Proof :  It  simply  follows  from  rank  V  =  r.  Hence  there  exist  r  linearly 
independent  rows  in  V.  From  AV  =  V  K  ,  where  f(K)  =  a  it  follows  that 

V  =  PV. 


I 
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This  result  says  that  with  the  mere  reordering  of  states,  and 
likewise  reordering  of  rows  of  V  there  can  always  be  found  inverible, 
i.e.  the  solution  of  the  Riccati  equation,  now  for  the  redefined  system 
(x,A) .  The  singularity  of  is  related  to  the  controllability  and 
observability  of  the  subsystems  of  A.  From  the  form  of  and  ^  (2-9) 
and  (2.10),  in  which  L  appears  once  as  a  controller  feedback  matrix  for 
the  pair  (Ai^A^)  and  the  other  time  as  an  observer  feedback  matrix  for 
the  pair  ,  it  can  be  seen  that  a  necessary  condition  for  to 

be  invertibel  (i.e.  for  L  to  exist)  is  that  the  modes  of  uncontrol¬ 
lable  through  A^2  belong  to  a =  o(B^),  and  the  modes  of  A22  unobservable 
through  A^0  belong  to  the  complementary  spectrum  a  =  cr^)-  This  becomes 
apparent  by  noting,  for  example  in  the  case  of  uncontrollability,  that  if 
a  ^  is  an  uncontrollable  mode  of  A^,  then  for  any  L  this  X  will  be 
in  c(B^),  which  is  contrary  to  the  requirement  that  X€ct^  =  cr^). 

It  should  be  noted  that  in  most  applications  the  approach  via  reordering 
rows  of  V  to  assure  invertibility  of  is  viable  one.  However,  in  the 
case  of  standard  Riccati  equation,  because  of  the  structure  of  E,  it  is 
not  possible  to  permute  rows  of  V,  except  within  and  V ^  in  which 
case  it  does  not  change  singularity  of  V^.  However,  under  usual  assumptions 
on  controllability-observability  it  is  known  that  the  unique  solution  for 
L  exists,  which  means  that  is  invertible.  For  some  algorithms  for  the 
time  scale  decomposition,  one  of  which  will  be  presented  in  the  next 
section,  it  is  important  to  make  reordering  of  states  to  assure  inverti¬ 
bility  of  before  the  algorithm  is  applied.  From  the  form  of  the 
transformation  (2.11)  and  the  requirement  that  X(B^)  <  X(B9)  it  can  be 
concluded  that  the  reordering  of  states  should  be  such  that  x^  contains 
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0 


the  slow  variables  of  the  model  and  x2  predominantly  fast  variables.  For 
cases  in  which  physical  considerations  or  previous  experience  cannot  help 
make  such  ordering .states  can  be  ordered  so  that  the  first  row  of  A  has  the 
smallest  norm,  the  second  row  has  the  next  largest  norm,  and  so  on. 


2.3.  Practical  Aspects  of  the  Spectrum  Separation  Problem 

In  practice  R(L)  cannot  be  made  zero,  except  in  very  special  cases. 
Therefore,  it  is  of  interest  to  see  how  accurately  ct(B^)  and  ct(b7)  of 


!  B 


A  = 
v 


1  12 
R(L)  B2 


(2.26) 


in  which  R(L)  is  supposedly  small,  reproduce  and  cr  .  Equation  (2.26) 
represents  the  lefthand  side  of  (2.7)  for  some  L;  we  assume  it  is  close 
to  the  solutuon  corresponding  to  a  given  c^.  The  answer  to  the  above 
question  is  given  in  the  following  theorem,  given  in  [32]. 

n  — 

Theorem  2.10:  Let  !|a!|  =  trace  (A*A)  and  define  5  =  inf  { i'PB^-B^Pii : 

P€R(n-r),<r .  P 


,  iIpiI  =  1}.  if 
|R(L)||  1!a 


12 


<}. 


r)xr 

then  there  exists  a  unique  P€R  satisfying 


lit 


P,i  <  2 

0 


(2.27) 


(2.28) 


such  that 


Im 


'^1 


(2.29) 


is  an  invariant  subspace  of  A  .  Moreover  ^(A^)  is  the  disjoint 


union 


ct(Av)  =  cr^+A^P)  U  a(B2  -  PA12)  . 


(2.30) 


In  practical  applications,  when  R(L)  is  sufficiently  small  we  approximate 
|"B1  A12 

A  with  ,  so  that  ct(B.)  approximates  o(B-  +A.„P)  and  cr(B_) 

V  0  B2  1  i  z 

approximates  cr^-PA^),  and,  finally,  Im  *  approximates  Im  *  .  The 

theorem  states  that  the  error  is  directly  proportional  to  P,  whose  norm 
is  bounded  by  the  norm  of  R(L).  In  other  words  there  is  a  continuity  in 
the  change  of  spectrum  and  the  invariant  subspace  for  small  changes  in 
R(L) .  The  number  6  in  the  theorem  is  a  measure  of  spectrum  separation 
between  B^  and  B2 .  If  B^  and  B0  are  symmetric  this  number  is  equal  [32] 
to  the  smallest  in  magnitude  difference  between  any  two  eigenvalues  of 
B^  and  B2 .  This  number  is  positive,  except  when  B^  and  B2  contain  common 
eigenvalues,  in  which  case  it  is  equal  to  zero.  However  under  Assumption  2.1 
it  will  always  be  positive.  Condition  (2.27)  indicates  that  the  Riccati 
equation  has  to  be  solved  more  precisely  if  either  there  is  a  strong 
connection  between  B^  and  B2  subsystems  through  A^2  or  spectrum  separation 
between  the  two  subsystems  is  small. 

We  now  use  the  theorem  to  give  an  alternative  interpretation  of 
some  of  the  standard  results  in  the  singular  perturbation  theory  [8]. 

When  the  state  matrix  of  (1.9),  which  is 
!  A11  A12  1 

A'  =  I  .  '  is  transformed  by  means  of  the  transformation  (2.11), 

I  21  22  ■ 

I_  t  e  ! 

A2?  A?i 

the  corresponding  Riccati  equation  is  R(L)  =  -  LA.^  +  ~jT~  L  -  LA^L  +  -g— 

For  small  £  the  dominant  terms  in  the  equation  are  those  containing  £. 


Hence  a  natural  choice  of  an  approximate  L  solving  the  equation  is  the 
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n 


I 


L 


one  which  annihilates  ^-dependent  terms.  This  leads  to  L  =  ~A22A21’  an<* 

B^  -  A^  -  ^12^22^21  s  Ag,  which  is  the  slow  state  matrix  as  given  by 
(1.11).  To  apply  the  Theorem  2.10  we  have  to  check  the  condition  (2.27). 

It  can  be  shown  that  for  nonsingular  A 22  and  sufficiently  small  £,  6  =  0(i), 
and  hence  there  exists  an  £  such  that  for  all  £c(o,£  ]  the  condition 
(2.27)  is  satisfied.  Therefore,  there  exists  a  P  =  0(6),  such  that  the 
spectrum  of  A’  is  given  as  a  disjoint  union 

c(A')  =  <?(AS+A12P)  U  cr(-|^  -  A22A21A12  -  PA12)  ,  (2.30) 


and  the  slow  eigenspace  of  A'  is  spanned  by 


:A22A21  +  P 


Because  of 


P  =  0(C)  this  result  exhibits  the  continuity  of  the  slow  eigenvalues 
and  the  slow  eigenspace  with  £€(0,£*] .  A  slight  difficulty  in  extending 
the  argument  to  hold  for  £  =  0  is  avoided  in  [10],  where  in  addition  the 
continuity  of  the  inverse  of  the  fast  eigenvalues  and  the  fast  eigenspace 

ic 

with  £€[0,£  ]  is  proved  as  well. 


2.4.  Relation  Between  Singular  Perturbation,  Time  Scale  Decomposition 
and  Modal  Methods 

In  the  previous  section  it  has  been  shown  that  the  time  scale 
decomposition  and  the  modal  method  for  reduced  order  modeling  are  the 
same  in  regard  to  the  slow  subsystem.  In  addition  time  scale  decomposi¬ 
tion  formulation  reveals  the  fact  that  with  the  same  expense  spent  for 
obtaining  the  slow  subsystem  state  matrix,  the  fast  state  matrix  can  be 
obtained  as  well.  For  some  engineering  applications  it  is  not  necessary 
to  find  exact  invariant  subspaces,  or  equivalently  the  exact  solution  of 
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the  Riccati  equation  [11].  Instead,  an  iterative  algorithm  for  the  time 
scale  decomposition  is  established,  whose  results  after  only  few 
iterations  are  then  used.  In  this  section  we  will  relate  so  defined 
time  scale  decomposition  with  the  singular  perturbation  and  the  modal 
methods . 

A  basis  for  comparison  of  the  methods  is  found  in  the  simultaneous 
iteration  method  for  computing  invariant  subspaces  of  general  (nonhermitian) 
matrices  (SI  for  short)  [32].  We  will  first  review  relevant  facts  about 
this  method  and  then  establish  the  relation. 

For  matrix  A  let  X^,  X^,  ••  be  eigenvalues  ordered  so  that 

\\\  *  U2I  *  ...  2  |XJ  (2.31) 

Let  denote  the  spectrum  consisting  of  the  r  largest  in  magnitude 
eigenvalues,  and  let  S r  be  the  corresponding  r-dimensional  invariant 
subspace.  Furthermore,  let  be  r-dimensional  subspace  with  nonzero 
projection  on  along  the  complementary  invariant  subspace.  Then,  if 
l^rl  >  I ^r+}J  >  A^Q^  tends  to  Sr  as  k  tends  to  infinity.  For 
r  =  1  this  corresponds  to  the  well  known  power  method  for  obtaining  the 
eigenvector  for  the  largest  eigenvalue  X^. 

Before  presenting  an  algorithm  which  uses  this  sequence  to 
find  a  basis  of  the  invariant  subspace  let  us  recall  the  result  about 
Schur-decomposition  of  a  matrix  [32,38]. 

Theorem  2.11:  Let  A€Rn*n.  Then  there  is  a  unitary  x€Cn  n  such  that 

S  =  X*AX  (2.32) 

is  upper  triangular.  The  matrix  X  may  be  chosen  so  that  the  diagonal 
elements  of  S  which  are  eigenvalues  of  A  are  in  descending  order  of 


absolute  value. 
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For  distinct  eigenvalues  corresponding  columns  of  X  are  unique  up  to 

multiplication  with  a  unit  length  number;  for  multiple  eigenvalues  of 

multiplicity  p,  the  corresponding  p  columns  in  X  span  unique  p-dimensional 

subspace.  If  X  denotes  first  r  columns  of  X  then  A*Im  X  = Im  X  ,  i.e. 

r  r  r 

these  columns  span  r-dimensional  invariant  subspace  corresponding  to  o  . 
Column  vectors  in  X  are  usually  called  Schur  vectors. 

Now  let  us  explain  one  step  of  SI  algorithm  for  finding  a  basis 

k 

of  A  Q  .  Define  rxr  matrix  B  as 
r 

B  =  V*AVk  (2.33) 


and  apply  the  decomposition  of  Theorem  2.3  to  B 
* 

‘  \  ~B 


<x£)  B  XB  =  ST 


(2.34) 


B  k 

Let  Qk  =  V^X^,  which  is  obviously  a  basis  for  A  Q. 


A  basis  for  the  next 


step  is  obtained  as 


<k+l 


-  AQ, 


(2.35) 


To  summarize,  the  algorithm  consists  of  orthonormalizing  an  arbitrary  basis 
matrix  (Q  ,  say)  to  get  V,  then  forming  B  (2.33),  decoupling  low 

K.  <  1 

order  matrix  B  according  to  (2.34)  and  using  unitary  matrix  to  perform 
the  next  step  (2.35). 

Convergence  behavior  of  this  algorithm  is  given  in  the  following 
lemma  due  to  [32]. 

Lemma  2 . 12 :  Let  S  be  equivalent  matrix  to  A  according  to  Theorem  2.11. 

Let  E^  denote  the  matrix  consisting  of  columns  1  through  2  of  I  and  E9 

denote  the  matrix  consisting  of  columns  2+1  through  r.  If  J\J  >  |  J  and 

|\  |  >  |  a.  I  then  there  are  matrices,  W,  ?Rn"r  with  orthonormal  columns 
1  r  r+1  k 


such  that 


divided  on  Wfc  =  (W^Wj),  W  CR®**,  V  such  that 

1.  In(Wk)  =  SkQ 

(2. 

2.  V1  =  Ex  +  0(|Xr+1/Xi|k) 

(2. 

3.  W2  =  E2  +  0(|Xr+1/\r|k) 

(2- 

The  lemma  shows  the  speed  with  which  individual  columns  in  Q  of  (2.35) 
tend  to  the  corresponding  Schur  vectors.  The  speed  is  linear  with 
the  factor  for  the  i-th  Schur  vector, defined  as 


e 

i 


9 


(2. 


and  the  convergence  of  the  algorithm  is  global. 

Practical  details  about  the  realization  of  the  algorithm  are 
contained  in  [32].  For  our  purpose  we  need  only  the  established  result 
about  global  convergence  of  the  algorithm  with  the  linear  speed  defined 
by  . 

We  also  rewrite  the  algorithm  (2.33)  -  (?..34)  in  the  form 


Using  this  form  of  the  algorithm  we  now  characterize  the  convergence 
properties  of  a  class  of  iterative  algorithms  for  solving  Riccati 
equation  represented  by  (1.19)  [13  3  -  We  rewrite  it  here  for  convenience 

Vi"  \  '  <A22'LkA12)'lR<Lk)  <2' 

Theorem  2.13:  Let  r  be  given  such  that  |\  |  >  with  the  eigen- 

r  r»lki 

values  of  A  ordered  according  to  (2.31).  Let  V,  =  j  |  be  an 
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nxr  matrix  with  nonzero  projection  of  column  vectors  on  the 
r-dimensional  slow  subspace  of  A.  If  exists,  denote  =  V2k^lk  ' 
Then,  if  Lq  =*  Pq,  the  sequences  P^  generated  by  (2.40)  and  generated 
by  (2.41)  are  identical.  Furthermore,  |X(B^(L^))|  <  |  X(B^  (Lj^) )  1  for 
sufficiently  large  k. 

Proof:  Writing  (2.40)  in  expanded  form,  substituting  P^V^  for  V^, 
eliminating  from  the  equations  and  regrouping  the  terms,  directly 
results  in  P^  =  L^.  Since  Im  tends  to  the  slow  subspace  of  A,  due 
to  the  Lemma  2.11,  B^(L^)  will  eventually  contain  the  r  slow  eigenvalues 
of  the  system. 

This  result  shows  that  the  Riccati  iterations  (2.41)  have 
global  convergence,  with  the  speed  determined  by  the  largest  6^,  i.e.  6^ 
which  corresponds  to  the  slowest  converging  Schur  vector  in  the  SI 
method.  Furthermore  the  spectrum  of  B^  =  A^  +Ai2Lk  ten<*s  to  CTS>  which 
means  that  Lk  tends  to  .  In  [15]  it  is  shown  that  is  the  only 
stable  equilibrium  solution  of  the  matrix  differential  equation 


L  =  R(L) 


(2.42) 


The  following  relation  between  singular  perturbation,  modal  methods  and 
the  iterative  time  scale  separation  is  immediate  consequence  of  the  above 
theorem. 

Corollary  2.14: 

Let  B^(k)  be  defined  as  B^(k)  =  A^  +  Ai2  ^k’  where  L^ 
is  defined  by  (2.41)  for  L^  =  0.  Then 

(i)  B^(l)  is  the  slow  state  matrix  obtained  from  singular 
perturbation  approach  (1.11) 
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(ii)  B^(k),  k  -»  ®  is  the  modal  method  slow  state  matrix,  and 
(iii)  B.(k),  1  ^  k  <  00  ,  is  the  slow  state  matrix  of  time 
scale  separation. 

Proof :  Substituting  Lq  =  0  into  (2.41)  gives  =  -  A22  A2^  and  B^(L^) 
of  the  form  given  in  (1.11).  This  proves  (i).  From  Theorem  2.12 

Im  ^  tends  to  the  invariant  slow  subspace  of  A  as  k  goes  to  infinity, 

L\J 

which  due  to  the  Remark  2.8  means  that  B^(L^)  tends  to  the  slow  subsystem 
matrix  obtained  from  modal  methods.  The  case  (iii)  is  the  case  in 
between  singular  perturbation  and  modal  methods,  which  offers  advantage 
over  the  exact  modal  method  when  low  accuracy  is  acceptable  but  £^_  is 
not  sufficiently  small  that  only  one  iteration  of  the  singular  perturbation 
approach  sufficies. 

This  corollary  ,  in  view  of  Lemma  2.11  and  Theorem  2.13  gives  a 

numerical  interpretation  for  why  the  singular  perturbation  approach  is 

successful  in  case  of  small  £;  it  is  simply  because  the  speed  of  conver- 

—  k 

gence  is  so  high  (proportional  to  £  )  that  only  one  iteration  is  needed 

R  (L) 

to  make  sufficiently  small.  The  corollary  also  suggests  one  meaningful 

interpretation  of  the  small  parameter  £,  as  being  the  eigenvalue  ratio 


£  = 


r41 


We  now  give  several  examples  to  illustrate  properties  of  the 
iterative  time  scale  decomposition. 

Example  2.15:  In  this  example  we  demonstrate  the  importance  of  the  state 
ordering  for  the  time  scale  decomposition  .  Given  (2.43)  and  r=2, 


X  = 


5  0  1 
Oil 
0  0  1 


(2.43) 


it  is  evident  that  the  state  matrix  has  two  Jordan  blocks;  one  of 

size  1  corresponding  to  the  simple  eigenvalue  5  and  the  other  of  size  2 

corresponding  to  the  eigenvalue  1.  From  the  triangular  structure  of 

this  matrix  it  is  evident  that  the  spectrum  separation  is  already  achieved 

for  any  choice  of  a .  However,  the  time  scale  separation  is  not  achieved 

since  the  large  eigenvalue  5  appears  in  the  "11"  block-  We  specify 

r  =  2  and  want  to  find  matrix  L, ,  which  is  a  part  of  the  transformation 

a 

(2-11)  so  that  cr(s^)  =  fl,l}.  The  Riccati  equation  has  the  form 

R(D  =  -  -i  +  <v2>Q  °)  •  <v2>(!)  <v2>  ■ 0 

from  which  two  solutions  can  be  found  easily:  =  (0  0)  and  =  (-4  0). 

For  the  first  solution  has  eigenvalues  5  and  1,  i.e.  nothing  is  changed 
from  (2.43).  Only  the  second  solution  results  in  the  desired  spectrum 
of  B^,  i.e.  Ld  =  L, .  Application  of  the  iterative  algorithm  (2.19)  to  this 
system  results  in  the  convergence  behavior  given  in  fig.  2.1a.  After 
initial  divergence  the  algorithm  achieves  asymptotic  behavior  predicted 
by  the  Lemma  2.12.  To  start  the  algorithm  we  have  perturbed  slightly 
since  for  the  usual  initial  condition  =  A22  ^1  the  algorithm  would 
stay  at  the  unstable  solution  L  =  0.  This  example  shows  several  interesting 
points.  First,  it  emphasizes  the  importance  of  state  ordering,  for  if 
the  state  ordering  was  such  that  A^  had  the  eigenvalues  close  to  the 
eigenvalues  of  the  desired  B^,  the  algorithm  would  have  started  from  the 
falling  side  of  the  curve  a,  in  fig.  2.1.  The  situation  is  illustrated 
by  the  curve  b,  which  corresponds  to  the  ordering  of  states  (x2>X2>x^) 
and  a  small  A2^.  This  is  also  one  of  the  reasons  for  the  success  of  singular 


perturbation  approach,  in  which  it  is  required  that  the  states  be  ordered 
according  to  their  speeds  (1.9).  This  example  shows  that  even  though  the 
eigenvalue  separation  exists,  only  after  the  states  are  properly  ordered 
one  step  approximation  of  the  singular  perturbation  approach  is  successful 
as  indicated  by  the  curve  b  as  opposed  to  the  curve  a.  Second  character¬ 
istic  of  the  example  is  that  it  contains  multiple  eigenvalues,  and 
still  achieves  the  asymptotic  behavior  predicted  by  Lemma  2.12.  In  other 
words  the  individual  generalized  eigenvectors  corresponding  to  multiple 
eigenvalues  may  not  be  unique  but  the  space  they  span  is  unique,  as 
discussed  in  connection  with  Theorem  2.11. 

Example  2 . 16 :  With  this  example  we  demonstrate  applicability  of  the 
techniques  of  this  chapter  to  the  solution  of  the  standard  Riccati  equation 
(2.14),  associated  with  the  optimization  on  infinite  interval.  The 
problem  is  to  find  a  control  u  which  will  minimize  criterion  (2.16)  for 
the  system  with 


A  = 


0 

0 


1 

-2 


E  = 


0  i 


t  ,  R  =  1,  Q  = 


T 

=  C  C 


0  1  I 


It  is  easy  to  check  that  (A,B)  is  controllable  pair,  and  (A,C)  is 


observable  pair.  In  fact  the  system  is  stable  since  eigenvalues  are 
0  and  -2.  The  objective  of  the  design  is  to  improve  its  performance. 
The  matrix  E  (2.17)  for  this  system  is 
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m 


whose  eigenvalues  are  (-2, -1,1, 2).  If  we  define  spectrum  separation 

-1  x 

problem  for  E  with  cr ^  ~  (-2,-1),  then  corresponding  B^  =  A  -  BR  B  L 
would  be  the  closed  loop  feedback  matrix  whose  spectrum  would  be  ant^ 

-1  x 

the  feedback  gain  K  from  u  =  Kx  would  be  K  =  -  R  B  L.  The  spectrum  of 
B^  indicates  improvement  in  performance  over  the  open  loop  system. 

Now,  we  want  to  use  globally  convergent  subspace  iteration  method  for 
computing  the  solution  L.  Since  the  eigenvalues  of  are  not  the 
largest  in  magnitude  eigenvalues  of  E,  we  perform  bilinear  transformation 


E  =  (E  -  col)  (E  +cpl) 


-1 


(2.44) 


tr» 


I* 


which  has  the  property  to  preserve  the  invariant  subspace  corresponding 
to  O'  and  has  eigenvalues 


ME) 


X(E)  -co 


ME)  +9 

so  that  the  n  eigenvalues  corresponding  to  Re  ME)  <  0  can  be  made  to  be 
the  largest  in  magnitude  eigenvalues  of  E  with  a  suitable  choice  of  co. 

A  co  that  sufficies  is 


(2.45) 


co  >  sup{Re  X.  (E),  i  =  1,2,  ...,  2n} 


(2.46) 


In  this  example  we  estimate  cp  =  5 .  We  now  use  subspace  iterations  (2.40) 
with  E  and  initial  guess  randomly  selected  as 


V0  = 


1 

.1 

.3 

.5 


-.2 

1 

.7  ; 

i 

-.8 


After  each  iteration  k  we  compute  i;R(L^)",  and  give  it  in  fig.  2.2  as 


« 


Log  i|R 


t 


a  function  of  the  number  of  iterations.  The  behavior  is  as  predicted 

by  the  Lemma  2.11.  After  11  iterations  elements  of  Lin  differ  from  the 

r  i  11 


corresponding  elements  of  the  exact  L 


6  2 
.2  1 


by  less  than  0.01. 


Convergence  is  faster  for  systems  whose  closed  loop  system  matrix  has 
larger  margin  of  stability  since  in  this  case  the  eigenvalue  separation 
between  the  eigenvalues  with  positive  real  part  and  the  eigenvalues  with 

\h-i® 

negative  real  part  is  bigger,  and  hence  £•  =  - — —  is  smaller.  This  example 

VE) 

demonstrates  the  possibility  of  using  the  advantage  of  having  globally  con¬ 
vergent  algorithm  for  solving  the  Riccati  equation,  for  either  initializing 
the  faster  converging  algorithms,  such  as  [33],  which  require  initial  guess 
close  to  the  final  solution,  or  directly  use  it  for,  say,  large  sparse 
system  matrices  A.  Because  of  the  simplicity  of  the  algorithm  it  can  also 
be  used  for  smaller  order  systems  when  limited  computer  capabilities  are 
available . 

Example  2.17:  This  example  illustrates  time  scale  decomposition  of 

one  machine  infinite  bus  power  system  model.  The  form  of  the  model  and 

the  linearized  state  matrix  [34]  are  given  in  Appendix  A.  The 

iterative  time  scale  decomposition  using  the  algorithm  (2.41)  is  exhibited 

in  fig.  2.3d,  for  the  case  r  =  2,  which  is  the  case  with  the  fastest 

convergence.  The  behavior  of  the  subspace  iterations  for  the  same  model 

and  different  r  is  given  in  fig.  2.3.  The  convergence  criterion  for 

B  *  B  12  r 

the  latter  is  derived  from  (2.34)  by  defining  E  =  (X^)  BX  -Sg  = (e  e  ...  e  ). 
This  example  will  be  completed  in  the  next  section  after  the  block- 
diagonalization  is  discussed,  by  showing  a  relation  between  the  number  of 
iterations  in  the  time  scale  decomposition  and  the  quality  of  the 
approximation. 
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2.5.  Block  Diagonalization 

After  the  spectrum  separation  is  completed,  the  system  state 
matrix  is  in  the  block  triangular  form  (2.12).  For  some  application 
it  is  useful  to  have  the  two  subsystems,  whose  dynamics  is  contained  in 
and  B2 ,  completely  decoupled.  In  geometrical  terms  it  amounts  to 
finding  a  basis  for  invariant  subspace  corresponding  to  the  spectrum  cr^, 
as  indicated  by  (2.4).  In  matrix  terms  we  can  apply  a  similarity 
transformation  to  (2.12),  similar  to  (2.11) 


I  h! 


0  I 


x  =  T  x 
v  H  v 


(2.47) 


which  does  not  change  the  physical  meaning  of  the  last  n  -  r  components 
of  xv-  When  this  transformation  is  applied  to  (2.12)  it  results  in 


where 


T  Bi  V  T-i_  Bi  P(H) 

H  H 

0  B_  0  B„ 


P(H)  =  HB2  -  BXH  +  A12 . 


(2.48) 


(2.49) 


The  solution  of  P(H)  =  0  is  a  function  of  the  solution  L  of  R(L)  =  0. 

It  is  well  known  that  equation  (2.49)  always  has  a  solution  when 
Assumption  2 . 1  is  satisfied.  In  case  when  a =  c(B^)  contains  r  slowest 
modes  of  A,  this  solution  can  be  found  as  the  limit  of  the  sequence 


Hk+1  '  (BA  -  A12>B2  • 


This  can  be  seen  by  writing  (2.50)  as 


hk+l  =  Ahk  +  b  =  (I  ©  B2)-[(Bi  ©  Dhk-  a12) 


(2.50) 


(2.51) 
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where  (g>  denotes  the  Kronecker  product.  The  eigenvalues  of  A  are 
X (A)  =  X(B1)/X(B2)  so  that  |\(A) |  < 1  and  hence  converges  to  a  unique  solution. 

Combining  transformations  (2.11)  and  (2.47)  we  get  the  complete 
transformation  between  x  and  (y,z)  as 

y  I-HL  H  ^X1  '  . 

=  2  j  s  T_ix  (2.52) 

z  -L  I  X 

where  L  and  H  satisfy  (2.13)  and  (2.49),  respectively. 

For  some  engineering  applications  where  time  scale  decomposition 
is  required,  only  a  few  iterations  of  (2.41)  and  (2.50)  for  L  and  H  are 
sufficient.  Figure  2.4  is  the  power  system  of  Example  2.17  where  the 
responses  of  the  system  based  on  subsystem  integration  of  (B^,y) 
and  (B2»z),  combined  through  (2.52)  into  x  are  compared  with  the  slow 
approximation  y.  The  matrices  B^,B2  and  T  contain  L  and  H  obtained  after 
two  iterations  of  (2.41)  and  (2.50). 

2.6.  Time  Scale  Separation  in  Weakly  Nonlinear  Systems 

Simulation  is  still  dominant  method  of  analysis  for  many 
systems,  especially  for  nonlinear  systems.  An  objective  of  the  analysis 
is  to  investigate  the  effect  of  parameter  changes  on  system  behavior, 
or  to  detect  system  instability.  The  emphasis  is  usually  on  the  speed 
of  simulation  rather  than  on  accuracy.  Yet,  the  speed  of  simulation 
is  constrained  by  the  wide  range  of  dynamic  phenomena  encountered  typically 
in  large  scale  systems.  In  other  words  some  variables  of  the  system 
respond  much  faster  than  the  others  (in  linear  systems  this  is  to  say 
that  there  is  a  spread  of  eigenvalues).  We  say  that  system  possesses 
multitime  scale  property.  This  property  is  indicated  by  writing  the 
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system  model  in  the  singularly  perturbed  form 


i,  -  VVV 

(2.53) 

ei2  '  WV 

(2.54) 

for  x  in  some  domain  d,  as  in  (1.9)  for  linear  systems.  Assuming 
that  the  fast  subsystem  is  stable  along  any  trajectory  initiated  in  D, 

[8], the  analysis  of  (2.53),  (2.54)  is  performed  in  two  time  scales 
avoiding  so  the  stiffness  of  the  model,  in  much  the  same  way  as  in 
the  analysis  of  the  linear  singularly  perturbed  systems.  The  problem 
arises  when  all  states  have  mixed  fast  and  slow  components,  i.e.  in 
the  mixed  state  case.  In  this  case  meaningful  decomposition  of  the  state 
vector  as  in  (2 .53) , (2 .54)  is  impossible. 

This  section  addresses  the  problem  of  the  simulation  of  two 
time  scale  mixed  state  systems  transformable  to  the  form 

k  =  Ax  +  f(x),  f(0)  =  0  (2.55) 

1* 

if  f(x)£C  and  has  small  Lipschitz  constant.  We  will  call  such 
systems  weakly  nonlinear  since  basic  characteristics  of  their  behavior 
are  determined  by  the  linear  part.  Hence,  two  time  scale  property 
of  (2.55)  implies  that  the  matrix  A  has  two  groups  of  eigenvalues  widely 
separated  in  magnitude.  If  T  is  the  block  diagonalizing  transformation 
(2.52)  for  the  matrix  A,  such  that  slow  eigenvalues  are  in  and  fast 
eigenvalues  are  in  (B^  and  B ^  are  defined  by  (2.9)  and  (2.10)), 
then  a  transformation  of  variables 


is  a  class  of  functions  with  continuous  first  derivative. 
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x  =  T 


j  •  [:]■[;: 


(2.56) 


can  be  used  for  the  nonlinear  model  as  well,  in  which  case  (2.55)  becomes 


y 

Bi 

o' 

y 

.  =  _ 

_  0 

B2- 

_2  . 

+  Wf(x) 


(2.57) 


where  W  is  the  inverse  of  T  given  explicitly  by  (2.52).  In  the  above 
model  y  and  z  represent  predominantly  slow  and  predominantly  fast 
variables  respectively,  which  are  decoupled  through  linear  part,  and 
weakly  coupled  through  the  nonlinear  part  (the  weak  coupling  is  the 
consequence  of  the  assumption  on  f  as  having  small  Lipschitz  constant). 
The  system  (2.57)  is  now  in  singularly  perturbed  form  and  the  standard 
singular  perturbation  approach  can  be  applied: 

(i)  solve  the  slow  subsystem 


y  "  Bjy  +  wxf(y2z) 


0  =  b2z  +  w2f(y 2z) 


(2.58) 

(2.59) 


for  y  and  z ,  and 
(ii)  solve  the  fast  subsystem 

v<tUJ 


z  =  B2z  +  B2z  + 


) 


(2.60) 


for  z.  Then  approximate  the  solution  x  by  x  =  T  |  ^  ;  . 

_  j 

Implicit  assumption  in  this  algorithm  is  that  the  fast  subsystem 
Jacobian 

rf  -  b2 + «2  a 


(2.61) 


evaluated  along  y  and  z  has  ReX(T^)  <  0  and  z(0)  is  in  the  appropriate 
domain  of  attraction  [8].  Under  this  assumption  z  can  be  expressed 
as  a  continuous,  unique  function  of  y,  and  hence  a  unique  solution  of 
(2.58)  and  (2.60)  exists  in  y  and  z.  Solving  this  equations  can  answer 
the  question  of  slow  subsystem  stability,  and  complete  system  behavior. 

If,  however,  we  assume  that  system  is  stable  and  want  to  observe 
its  slow  behavior  only»for  different  initial  conditions  in  the  stability 
domain,  or  after  a  change  in  parameters,  then  a  useful  combination  of 
the  linear  and  nonlinear  simulation  can  be  achieved  if  the  system  is 
presented  in  the  form  (2.57).  In  this  case  after  some  time  tc>  such 
that  f(x)  is  sufficiently  small,  the  simulation  can  be  continued  on  the 
low  order  linear  subsystem  y  =  B^y.  We  illustrate  now  some  of  the 
potentials  of  this  approach,  and  some  differences  from  other  related 
approaches  to  simulation,  on  the  one  machine  infinite  bus  power  system, 
whose  model  is  given  in  Appendix  A.  Let  us  denote  that  model  as  x  =  cp(x). 
For  integration  we  use  the  predictor-corrector  method  of  [ 35 3  ,  but  do 
not-  investigate  which  numerical  method  is  the  best.  We  will  show  a  set 
of  figures  with  responses  of  x^,  which  is  typical  of  what  we  want  to  show. 

On  each  figure  there  are  two  responses:  solid  line  is  always  the  full 
system  response,  i.e  the  solution  of  x  =  q>(x).  The  dashed  line  corresponds 
to  different  approximations  of  the  slow  response,  each  of  which  will  now 
be  explained. 

First  the  system  is  linearized  around  the  stable  equilibrium, 

and  the  linearized  state  matrix  A  is  given  in  Appendix  A.  Then 

f(x)  =  'o(x)  -Ax.  For  the  block-diagonalization  of  A  we  specify  r=2, 

j  =  O'  and  use  T  of  (2.52),  in  which  L,  H,  B,  and  B~  are  as  obtained 
2  s  1  2 


after  two  iterations  of  the  algorithms  (2.41)  and  (2.50),  respectively. 

Now  the  model  is  in  the  form  (2.57).  Figure  2.5  gives  x^  as  obtained  by 
solving  (2,53)  and  (2.54)  with  e  set  to  zero.  The  error  is  unacceptable. 
Figure  2.6  shows  the  nonlinearity  present  in  the  slow  system  response, 
since  the  dashed  curve  is  the  response  of  y  =  B^y .  The  true  nonlinear 
slow  response  should  be  a  line  which  would  pass  in  between  the  weggles  of 
the  solid  curve.  Figure  2.7  gives  the  response  of  the  slow 
component  of  x^,  which  result  from  solving  (2.58)  and  (2.59)  for  y  and 
using  the  back  transformation  x^  =  y  (which  comes  from  x  =  T  [^j  ). 

The  improvement  of  this  approach  (application  of  singular  perturbation  to 
the  transformed  system)  over  the  application  of  singular  perturbation  to 
the  original  model  is  significant.  On  the  next  figure,  fig.  2.8 
the  slow  response  of  is  given  as  computed  by  using  nonlinear  model 
(2.58),  (2.59)  for  the  first  two  seconds  only,  and  then  for  the  rest  of 
the  time  using  linearized,  low  order  model  y  =  B^y.  In  addition  larger 
step  size  is  used.  Here  the  nonlinear  model  is  used  to  initiate  properly 
the  slow  linear  subsystem.  The  effect  of  nonlinearities  is  evident  from 
comparing  these  responses  with  fig.  2.6. 

In  some  resent  approaches  to  numerical  integration  of  singularly 
perturbed  and  descriptor  variable  system  [36]  it  has  been  proposed  to 
integrate  the  original  system  (2.53),  (2.54)  but  adjust  the  step  size  of 
integration  according  to  the  slow  component  of  the  local  error  (see  [50]  for 
numerical  integration  details) .  For  "filtering"  the  errror  in  order  to 
get  its  slow  component  a  matrix  very  much  related  to  W^(2.53)  is  computed, 
whenever  the  step  size  has  to  be  changed.  Based  on  the  power  system 
example  only,  it  seems  that  that  approach  may  give  more  accurate  response 


Application  of  the  singular  perturbation  to  the  nonlinear  model  of  one 
machine  infinite  bus  system:  exact  systems  response  vs.  its  slow 
approximation . 
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of  the  slow  dynamics  if  required,  then  it  is  possible  by  solving  (2.58) 
(2.59).  However,  slow  dynamics  of  (2.58),  (2.59)  permits  much  larger 
step  size  than  in  the  method  of  [ 36] •  This  is  useful  when  a  quick 
approximate  solution  is  satisfactory.  The  difference  in  accuracy  between 
the  two  methods  is  illustrated  for  the  case  when  both  are  forced  to  use 


very  large  step  size,  fig-  2.9. 


3.  PARTIAL  COHERENCY 


3.1.  Introduction 

Application  of  the  modal  method  for  the  reduction  of  the  multi¬ 
machine  electromechanical  model  of  power  systems,  has  been  tried  in  [16,18]. 
However,  more  successful  has  been  the  use  of  coherency.  Two  generators, 
with  responses  denoted  by  x^(t)  and  x_.  (t)  are  called  coherent  if 

x^(t)-x.(t)  =  const.  (3.1) 

Groups  of  machines  in  which  any  two  generators  satisfy  this  criterion  are 
called  coherent  groups  or  coherent  areas.  Once  coherent  areas  are  known  the 
reduction  of  the  model  order  is  achieved  by  replacing  all  the  coherent 
generators  by  one.  A  critical  step  however,  is  the  grouping  of  the  machines 
into  areas.  Coherent  machines  are  identified  either  from  actual  or 
simulated  machine  responses  [17,18],  or  by  an  algebraic  evaluation  of  the 
modes  present  in  the  linearized  model  [16,20,21,23].  Most  analytical  methods 
require  that  machines  be  coherent  throughout  the  duration  of  their  transients. 
We  propose  here  a  less  demanding  definition  of  partial  coherency.  It  may 
be  interpreted  as  a  requirement  that  the  equivalent  machines  of  the  areas 
represent  as  closely  as  possible  a  preselected  group  of  modes.  In  the  next 
chapter  it  will  be  shown  that  in  the  particular  case,  when  these  modes  are 
the  slowest  modes  of  the  system,  the  resulting  area  decomposition  is 
relatively  independent  of  fault  locations  and  loading  conditions.  In  our 
approach  partial  coherency  is  related  to  the  spectrum  separation  of  the 
previous  chapter.  In  the  ideal  partial  coherency  case  the  solution  L  of  the 
suitably  formulated  Riccati  equation  has  the  elements  made  of  zeros  and  ones. 


which  unambiguously  define  coherent  areas.  The  subsystem  reproducing  comple¬ 
mentary  subspectrum  turns  out  to  model  intermachine  dynamics,  showing  that 
in  the  case  of  ideal  partial  coherency  the  interarea  and  intermachine 
dynamics  are  decoupled.  In  a  nonideal  case  our  approach  is  to  search  for  a 

solution  L  which  can  be  approximated  by  the  matrix  of  zeros  and  ones.  In 

this  case  areas  are  made  of  machines  which  are  near-coherent  in  the  pre¬ 
selected  modes,  and  weakly  coupled  rather  than  decoupled  area  and  intermachine 
dynamics . 

In  connection  with  power  system  stability  study,  there  is 
continuing  interest  in  the  development  of  a  systematic  area  decomposition 
procedure.  Our  grouping  algorithm  reduces  the  decomposition  procedure  to  the 
computation  of  a  basis  for  the  invariant  subspace  corresponding  to  the  given 
subspectrum,  and  a  Gaussian  elimination  of  a  low  order  matrix. 

In  the  next  section  we  present  an  electromechanical  model  and 
review  some  of  its  properties.  An  overview  of  existing  approaches  to 

coherency  will  be  given  in  Section  4.3.  In  Section  4.4  the  concept  of  partial 

coherency  will  be  introduced  and  its  properties  analyzed  for  ideally  decom¬ 
posable  systems.  A  grouping  algorithm  developed  for  near  decomposable 
systems  will  be  given  in  Section  4.5.  The  algorithm  will  be  illustrated  by 
a  48  machine  example.  Physical  meaning  of  the  variables  associated  with 
the  decomposition  of  a  system  into  coherent  areas  will  be  given  in  Section  4.6. 
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3.2.  Electromechanical  Model 

The  electromechanical  model  of  a  power  system  with  n  generators 
labeled  from  1  to  n  and  m  load  buses  labeled  from  n+1  to  n+m  is  given  by  [52] 

5 .  =  u>. 

l  l 


0 


M.oi.  +  D.oj  =  Pm  -  P  . 
i  x  i  i  i  gx 

N  , 

P„.  =  Z  v.  v.B. .sin(5.-6 . )  +  v  G 
gx  j=i  i  ]  ij  i  j  i  ii 

j#i 


i=l,2 . n 


N  2 

=  Z  v. v.B . .sin(6 . -6 . )  +  v.G. .  i=n+l,...,N 

Hi  j=1  i  j  i]  l  j  l  n 

j^i 


(3.2) 


(3.3) 


where 


N  =  n+m 

<5  ^  =  for  l<i<n  rotor  angle  of  machine  i;  for  n+1  <  i  <  N  bus  angle  of 
the  load  bus  i-n  (radians) 

=  speed  of  machine  i  (per  unit) 

Pnu  =  mechanical  input  power  (per  unit) 

Pg_^  =  generated  electrical  power  of  machine  i  (per  unit) 

=  moment  of  inertia  of  machine  i 
=  Damping  constant  of  machine  i  (per  unit) 

P^  =  negative  of  load  consumed  at  the  load  bus  i-n. 

B^  =  reactive  part  of  the  admittance  connecting  buses  i  and  j 
G  =  i-th  diagonal  entry  of  the  real  port  of  admittance  matrix. 
Assumptions  associated  with  the  model  are: 

Assumption  3.1:  Mechanical  input  power  Pm^  and  load  P  are  constant. 
Assumption  3.2:  Active  power  losses  are  neglected,  i.e.  active  part  of 
admittance  matrix  is  neglected. 
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Assumption  3.3:  Active  and  reactive  power  flow  are  decoupled,  i.e.  bus 

voltages  are  assumed  constant  in  the  model  (3.2),  (3.3). 
Model  (3.2)  and  (3.3)  contain  both  dynamic  and  static  equations.  In  linear 
system  theory  such  models  are  known  as  descriptor  variable  systems  [9,10]. 
Study  of  the  descriptor  variable  system  obtained  from  (3.2)  and  (3.3)  by 
linearization  around  a  stable  equilibrium  point  reveals  to  a  great  extent  the 
intermachine  oscillation  behavior.  The  linearized  model  is 


where 
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(3.5) 


(3.6) 


B..  =  reactive  part  of  the  admittance  between  buses  i  and  j. 
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If  exists  (by  the  assumption  on  stability  of  the  equilibrium, is  in 
fact  positive  definite)  (3.3)  can  be  reduced  to 


0  I 

-M-1K  -M_1D 


(3.7) 


where 


K  =  H  -H  „H„X  =  (k..) 
gg  gSL  U  Ig  ij' 


(3.8) 


Matrix  K  has  several  useful  properties  which  can  be  deduced  from  the  pro¬ 
perties  of  the  Jacobian  H.  First,  due  to  (3.6)  H  and  K  are  symmetric.  From 

(3.6)  it  can  be  seen  that  matrix  H  has  one  zero  eigenvalue  with  the  corre- 

T 

sponding  eigenvector  v  =  (1 ,1 , . . . , 1) .  The  same  is  true  for  the  matrix  K, 
namely  Kv^  =  0,  where  v^  is  the  n-vector  containing  first  n  components  of  v. 
This  property  can  be  derived  as  follows.  Let  v^  be  n-vector  containing  last 


n  components  of  v.  Then,  from  Hv = 0  we  have 


V2  ' 


(3.9) 


so  that 


Kv.  =  (H  -H  XX  ) v  =  H  v,  +  H  .v,  =  (H  H  .)v  =  0.  (3.10) 

1  gg  gl  U  ig'  1  gg  1  g«  2  gg  gl 

This  property  shows  that  each  diagonal  element  of  K  equals  the  negative  sum 
of  the  off-diagonal  elements  in  the  same  row.  Under  the  assumption 
Assumption  3.4;  [  5^— 5^  |  <  tt / 2,  if  we  have  that  the  off-diagonal 

elements  of  H  are  negative,  which  due  to  the  zero  eigenvalue-eigenvector 

property  means  that  the  diagonal  elements  are  positive.  The  same  is  true  for 

“1  * 

K  since  then  has  all  its  elements  positive  so  that 


Matrix  H,,  is  Minkowski  matrix  for  which  Theorem  A. 2  of  [57] 
applies  asserting  that  the  elements  of  are  positive. 
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is  negative,  and  hence  all  the  off-diagonal  elements  of  K  are  negative  as 
well.  Due  to  (3.10)  it  means  that  the  diagonal  elements  are  positive. 

Based  on  the  properties  of  matrix  K=  (k^)  it  can  be  assumed  that 
the  elements  k„  are  of  the  form  (3.6).  In  this  way  we  define  an  equivalent 
admittance  matrix  connecting  generator  buses  only.  Under  Assumption  3.4  the 
off-diagonal  elements  of  the  reduced  admittance  matrix  are  of  the  same  sign, 
and  its  diagonal  elements  are  negative  sums  of  the  off-diagonal  elements  in 
the  same  row.  An  alternative  way  to  get  the  same  kind  of  model  is  to 
transform  each  load  to  an  equivalent  admittance,  and  then  eliminate  all  load 
buses.  In  summary  we  will  assume  that  the  electromechanical  model  is  in  the 
form  (3.7)  in  which  k_  elements  have  the  form  (3.6)  for  some  equivalent 
admittance  elements  B  „  .  This  model  has  several  properties  which  allow  even 
further  simplifications. 

At  6*  and  m*,  the  eigenvalues  of  (3.7)  are  of  the  following  three 
types  (see  Figure  3.1  for  an  example): 

1.  a  zero  eigenvalue  corresponding  to  the  motion  of  all  the  machine 
angles , 

2.  a  small  negative  real  eigenvalue  corresponding  to  the  aggregate 
speed  of  all  the  machines,  and 

3.  (n-1)  pairs  of  lightly  damped  oscillatory  modes  which  typically 
range  in  frequency  from  */2  to  2  Hz. 

Models  involving  more  details,  such  as  excitation  systems  and  governors  would 
still  contain  the  above  set  of  eigenvalues,  modified  mostly  in  the  damping  and 
not  in  frequency  [19].  Since  the  small  damping  constant  does  not  signi¬ 
ficantly  affect  the  frequencies  of  the  oscillatory  modes  [27],  it  can  be 
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neglected.  See  Figure  3.1  for  a  typical  pattern  of  eigenvalues  of  (3.7). 

Thus  the  model  used  in  our  approach  is 

x  =  -M^Kx  =  Ax.  (3.11) 

The  properties  of  the  A  matrix  are  as  follows : 

(PI)  A  has  zero  eigenvalue  whose  eigenvector  is 

v*  =  (1  1  ...  1).  (3.12) 

Property  (PI)  follows  from  Av^  =  0,  which  is  due  to  (3.6)  as  the  sum  of  each 
row  in  A  ■  (a  )  is 

n  1  n 

I  =  M  Z  k . .  =  0 ,  i = 1,2, . . . ,n.  (3.13) 

j=l  ^  1  j-1  *3 

(P2)  A  is  diagonalizable  because  it  is  similar  to  the  symmetric 

matrix 

-M",/2KM",/2  (3.14) 

i/ 

where  M  2  is  the  positive  symmetric  square  root  of  M.  Thus,  all  the  eigen¬ 
values  A^  of  A  are  real.  Furthermore,  under  assumption  that  the  equilibrium 
of  (3.2)  and  (3.3)  is  stable,  they  are  all  nonpositive  [51]. 

(P3)  If  A  is  the  state  matrix  of  the  system  (3.7),  with  damping 
neglected,  then  the  relation  between  eigenvalues  and  eigenvectors  of  A  and 
A  is 

A  =  +  /Aj\  (3.15) 

If  x  is  an  eigenvector  of  A  corresponding  to  A^,  then  the  two  eigenvectors  of 
A  are  (x^, +  /a  xV.  Because  of  such  simple  relations  between  eigenvalues 
and  eigenvectors  of  2nx2n  matrix  A,  and  nxn  matrix  A,  we  choose  to  deal  in 
our  approach  only  with  the  lower  order  A  matrix. 
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There  are  two  ways  of  representing  disturbances  in  such  linearized 
models.  In  [19,21,23]  the  disturbance  is  modeled  as  an  external  input  to 
either  unreduced  model  (3.4)  or  the  reduced  (3.7).  All  the  standard 
disturbances  -  generator  outage,  line  switching  and  load  shedding  -  are 
modeled  as  the  equivalent  change  in  load  and  generation.  The  duration  of  the 
external  input  equals  the  duration  of  the  fault.  Hence  the  model  is 


x  =  Ax  +  Bu 

x  =  Ax 


0  <;  t  £  t 


t  >  t 


where 


A  is  the  state  matrix  of  (3.7)  and  B  = 


H 


g8- 

I 


«u>'1 


u  =  input ;  u  = 


AP 

m 

lap* 


(3.16) 

(3.17) 


=  "22"  block  of  the  Jacobian  H  modified  in  structure  to  account  for  a 

line  outage;  otherwise  equal  to  H^. 

t£  *  clearing  time  -  duration  of  the  disturbance. 

In  our  approach  disturbance  is  modeled  as  an  initial  condition  at  t=t  for 

c 

the  linearized  system  (3.11)  around  a  postfault  equilibrium  point.  In  other 
words  only  the  low  order  version  of  (3.4)  is  considered. 


3.3.  Review  of  Existing  Coherency  Methods 

Coherency  has  been  a  fruitful  approach  to  the  model  order  reduction 
of  power  systems.  The  need  for  repeated  simulation  of  larger  and  larger  power 
systems  has  emphasized  even  more  the  role  of  this  technique.  A  reflection  of 
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this  situation  is  the  existence  of  a  number  of  different  approaches  to 
coherency  identification.  In  this  section  we  give  an  overview  of  several 
characteristic  approaches,  in  order  to  identify  the  limitations  and  gain 
insight  into  the  problem. 

In  [37]  the  conditions  for  coherency  of  a  group  of  machines  are 
derived  starting  from  the  definition  (3.1)  and  using  (3.2),  (3.3),  and,  in 
addition,  the  equations  for  the  reactive  power  flow.  The  conditions  are 
expressed  in  terms  of  machine  internal  voltages  and  interconnecting  admit¬ 
tances  as  follows.  A  group  of  r  machines  in  an  n  machine  system  is  coherent 
if  and  only  if 


B.  .  V 

a.  -77^-  (V.  cos  6  -jV.  sin  6  )  =  ~  B  . 

M.  1  ir  J  1  ir  M  rj 
1  r 


i  =  l,2,...,r;  j  =  r+l,...,n  (3.18) 

b.  the  disturbance  is  outside  the  coherent  area. 

This  condition  indicates  how  strict  requirements  are  to  be  satisfied  for 
coherency  to  exist.  One  case  when  they  are  satisfied  is  when  there  are  several 
machines  on  the  same  bus.  To  deal  with  the  practical  systems  a  less  restric¬ 
tive  condition  is  needed.  The  one  of  reference  [37]  is  given  using  the 
linearized  model  (3.16).  If  r  machines  are  to  be  coherent,  and  if  the 
corresponding  states  are  the  first  r  states  of  x  in  (3.16)  then  the  condition 
for  coherency  is  expressed  as 


a.  II A  II  <  e  (3.19) 

b.  II  B^ul!  <  y 


for  some  small  numbers  £  and  y.  Matrices  and  are  rx(n-r)  ,,12"  block 
of  A  and  r*l  block  of  B.  In  other  words  coherency  is  a  result  of  weak 
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coupling  of  the  group  of  r  machines  with  the  rest  of  the  system  and  the 
disturbance  is  outside  the  group.  The  reference  [37]  does  not  contain  an 
algorithm  for  identifying  areas. 

In  [20,21]  coherency  is  related  to  controllability  observability 
of  (3.16).  The  output  is  defined  as  follows.  If  a  coherent-  area  is  to 
contain  r  machines,  that  is 


y 


Tr-l,r-l 


-1 

-1 


0 


x  =  Cx 


(3.20) 


where  x  has  been  ordered  so  that  first  r  components  are  the  states  of  the 
coherent  machiens ,  then  a  necessary  and  sufficient  condition  for  coherency  is 
that 

In  Q  CKer  C  (3.21) 

c 

where  0^  is  the  controllability  matrix  (B,AB,...,An  '*'B)  .  This  condition 
simply  says  that  from  any  of  the  machine  differences  in  the  coherent  area 
the  disturbance  should  not  be  seen,  either  because  the  controllable  modes  are 
unobservable  through  C  or  because  the  observable  modes  are  uncontrollable 
through  B.  In  the  case  of  nonideal  coherency  it  is  suggested  to  measure  the 
distance  between  the  two  subspaces  of  (3.21)  and  define  coherency  when  the 
distance  is  small.  An  algorithm  for  checking  the  condition  (3.21)  would  be 
based  on  singular  value  decomposition  [44]  of  several  matrices  of  the  same 
order  as  the  order  of  the  system,  and  is  therefore  inappropriate  for  large 
power  systems. 

The  approach  of  [23]  implicitly  exploits  the  controllability- 
coherency  relation.  Two  basic  ideas  are  realized.  First,  that  machines  which 
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are  coherent  in  the  model  (3.16),  i.e.  during  the  fault,  will  be  coherent 
throughout  the  entire  response  period.  Second,  coherency  is  essentially 
defined  with  respect  to  the  modes  which  are  excited  by  the  disturbance. 
Namely,  eigenvectors  of  the  modal  matrix  M  from 


x  ( t )  =  Mz  (t) 


(3.22) 


are  multiplied  by  the  initial  values  of  mode  excitation,  i.e.  the  matrix 


M  =  Mz  (0) 


(3.23) 


where  z(0)  is  2nx2n  diagonal  matrix  with  the  ith  diagonal  entry  equal  to  z.(0) 
is  considered.  Then  two  machines  i  and  j  are  declared  coherent  if  the  dif¬ 
ference  || nu -null  <  £  for  some  small  e,  where  m^  and  nu  are  row  vectors  of  M. 

The  method  requires  computation  of  all  eigenvalues  and  both  right  and  left 
eigenvectors  of  the  matrix  A.  For  large  power  systems  this  may  require  con¬ 
siderable  computer  time  and  memory. 

Instead  of  normalizing  all  the  eigenvectors,  as  in  (3.23),  the  idea 
of  using  only  a  few  preselected  modes  (a  few  eigenvectors  of  M)  to  define 
coherency  has  been  explored  in  [26].  There  the  coherency  of  machines  i  and 
j  with  respect  to  mode  k  is  defined  as  the  requirement  that 


angle (nu^nu^)  <  e  <  tt/2 


(3.24) 


where  m  and  m.,  are  (in  general)  complex  elements  of  the  matrix  M.  For  a 

IK  J  K 

real  eigenvector  this  requirement  means  that  the  elements  of  the  k-th  eigen¬ 
vector  corresponding  to  the  i-th  and  j-th  machine  have  to  be  of  the  same  sign. 
The  criterion  means  that  the  machines  i  and  j  are  accelerated  (decelerated) 
by  the  power  from  the  rest  of  the  system  for  more  than  half  a  period  of  the 
k-th  frequency.  Coherency  with  respect  to  the  group  of  modes  is  defined  so 
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that  the  criterion  (3.24)  is  satisfied  with  respect  to  each  mode  of  the 
group . 

The  idea  of  using  coherency  for  model  order  reduction  has  found 
overwhelming  response  from  industry  after  it  had  been  promoted  in  [41],  and 
in  particular  after  a  complete  computer  program  has  recently  been  developed 
and  its  utility  demonstrated  on  a  number  of  large  systems  [19].  The  method 
employed  is  a  direct  application  j f  the  definition  (3.1),  which  means  that 
to  obtain  system  response,  the  integration  is  needed  for  each  disturbance  and 
then  comparison  of  responses  for  each  t  on  some  time  interval.  A  fast  inte¬ 
gration  method  (trapezoidal  integration)  is  used  to  keep  the  numerical  burden 
on  the  reasonable  level.  The  method,  however,  does  not  provide  any  insight 
into  the  relation  of  coherent  areas  to  the  system  parameters,  as  in  the 
previous  methods. 

In  all  the  previous  methods  coherent  areas  depend  on  the  location 
of  the  fault.  Therefore  for  each  new  fault  reevaluation  of  the  areas  has  to 
be  made.  In  [39]  an  attempt  has  been  made  to  overcome  this  difficulty,  and  to 
come  up  with  the  areas  which  can  be  used  for  more  disturbances.  A  solution  is 
found  in  a  probabilistic  approach.  Namely,  a  set  of  disturbances,  i.e.,  inputs 
u  in  the  model  (3.16)  with  given  probabilities  is  used.  A  measure  of 
coherency  between  machines  i  and  j  is  computed  as 


JT(x . -x . ) 2dt  = 


TP  o  '  1 


Se .  . 
1  j  iJ 


where  p  is  selected  to  make  c..  finite  for  T-*-00, 

ij 


1  1  T 

s  =  —  /  E ; x ( t )  x(t) }dt. 


TP  o 


(3.25) 


» 


(3.26) 
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and  e„  is  the  vec-tor  of  zero  elements  except  i-th,  which  is  1  and  j-th 
which  is  -1.  It  is  easy  to  see  that  smaller  values  of  c^  indicate  that 
machines  i  and  j  are  coherent.  Of  particular  interest  is  the  case  of  the 
infinite  time  interval.  Then  S  can  be  expressed  as^ 

S  -  ATB(R+  pyT)BT(AI)T  (3.27) 

where  A  and  B  are  defined  by  (3.16),  u  is  the  expected  value  of  u,  and  R  is 
its  covariance  matrix.  A*  is  defined  as 


A 


I 


(3.28) 


where  is  a  2n-2  diagonal  matrix  of  nonzero  eigenvalues  of  (3.16)  and  X 

is  its  matrix  of  right  eigenvectors.  For  zero  mean,  independent  identically 

T 

distributed  disturbance  R+yu  =  I  and  S  of  (3.27)  depends  only  on  the  para¬ 
meters  of  the  system.  This  dependence  is  made  explicit  in  [39]  by  writing 
part  of  (3.27)  in  terms  of  transformed  matrix  A,  rather  than  via  eigenvalues 
and  eigenvectors.  However,  by  using  the  above  form  to  eliminate  all  the 
modes  of  the  system  which  cause  insignificant  changes  in  S,  it  has  been 
concluded  in  [39]  that  the  eigenvalues  of  the  reduced  system,  based  on 
coherency  and  eigenvalues  of  the  reduced  system  based  on  the  modal  method 
[15],  are  approximately  the  same. 

This  important  connection  between  coherency  and  modal  methods  has 
also  been  indicated  in  [16]  under  the  special  condition  that  L  of  (2.11)  has 


This  expression  is  different  from  [39]  in  that  it  is  independent 
of  the  choice  of  reference  machine  whose  choice  may  be  a  problem,  but  is  in 
fact  irrelevant  for  the  existence  of  coherent  areas  as  will  be  shown  later  in 
the  chapter. 
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the  structure  in  which  each  row  has  only  one  zero  element  equal  to  1,  and  all 
other  elements  equal  to  zero. 

Most  of  the  presented  coherency  methods  concentrate  on  the  develop¬ 
ment  of  an  alternative  coherency  criterion,  in  order  to  avoid  direct 
simulation,  implicitly  required  in  Definition  3.1.  Except  for  [19]  none  of 
the  methods  offer  an  algorithm  for  identifying  coherent  areas,  which  is 
efficient  enough  to  be  applied  to  large  power  systems.  Most  of  them,  except 
[39,26]  would  give  areas  which  are  valid  for  only  one  disturbance. 

The  coherency  definition  and  the  corresponding  algorithm  for 
identifying  coherent  areas,  which  will  be  given  in  the  rest  of  this  chapter, 
are  aimed  at  overcoming  the  difficulties  of  the  existing  methods  and 
retaining  their  good  characteristics.  We  use  the  result  of  [16]  directly 
to  define  the  so-called  ideal  coherency.  Systems  with  ideally  coherent  areas 
serve  to  give  us  an  insight  which  we  then  use  in  analysis  of  realistic  power 
systems,  in  much  the  same  way  as  the  descriptor  variable  system  analysis  of 
(1.11)  serves  to  give  insight  into  the  behavior  of  singularly  perturbed 
systems  (1.9).  A  result  of  this  approach  is  the  new  algorithm  for  identifying 
coherent  areas,  designed  for  application  to  large  power  systems.  Analytical 
study  reveals  that  so  obtained  coherent  areas  are  almost  independent  on  fault 
location  and  nonlinearities  of  the  system. 


3.4.  Partial  Coherency 

Considering  the  definition  of  coherency  (3.1)  and  the  model  (3.17) 
it  is  clear  that  coherency  can  be  achieved  only  if  the  initial  condition 
x(0)  is  such  that  only  r  modes  are  excited,  for  some  r<n.  This  is  also 
clear  from  the  coherency  condition  (3.21),  for  if  dim!m(Qc)  =  n,  this  would 
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mean  that  there  exists  a  single  disturbance  which  would  excite  all  the  modes, 
but  then  C  must  be  zero  matrix  which  means  that  there  is  no  coherency.  In 
other  words,  each  case  of  coherency  is  a  partial  coherency.  In  this  section 
we  will  give  a  new  formulation  for  the  existence  of  partial  coherency  in  terms 
of  generalized  Riccati  equation.  The  coherency  will  be  connected  with  the 
spectrum  separation  problem  treated  in  Section  2.2.  Let  us  first  define 
partial  coherency. 

Definition  3.5:  Given  r  eigenvalues  of  A  in  (3.11),  a  .  Then  machines  ”i" 
and  "j"  are  coherent  with  respect  to  cr^  (partially  coherent)  if  for  all  t  of 
interest,  possible  t£  [O,00),  and  any  initial  condition  x(0),  their  angles 
x^t)  and  Xj  (t)  satisfy 

xi(t)-xj.  (t)  =  z^(t)  (3.29) 

where  z..(t)  contains  none  of  the  modes  from  a  .  A  coherent  area  consists  of 
ij  r 

all  the  machines  coherent  to  each  other. 

We  note  that  in  this  definition  no  machines  from  different  areas  can 
be  coherent,  that  is,  no  coherent  area  can  be  divided  into  more  areas. 

Although  Definition  3.1  does  not  require  that  the  number  of  coherent 
areas  equal  the  number  of  modes  r,  systems  with  this  property,  which  will 
be  called  r-decomposable  systems,  are  of  particular  importance.  The  study  of 
r-decomposable  systems  is  an  essential  step  toward  the  analysis  of  more 
common  "near-decomposable"  systems,  that  is  systems  with  near-coherent  rather 
than  coherent  areas. 

Definition  3.6:  The  machines  "i"  and  "j"  are  near-coherent  if  in  Definition 

3.1  the  contribution  of  the  modes  from  o  in  z..(t)  is  small.  A  near-coherent 

r  ij 
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area  consists  of  all  machines  which  are  near  coherent  to  each  other.  An  r- 

near-decomposable  system  consists  of  r  near-coherent  areas. 

Our  approach  to  the  identification  of  coherent  areas  is  to  first 

consider  this  problem  for  r-decomposable  systems.  We  show  that  in  this 

idealized  case  the  solution  of  the  matrix  Riccati  equation  (2.13),  which 

2 

separates  the  spectrum  of  A  into  a and  o^,  automatically  groups  the  machines 
into  areas.  We  then  use  this  result  to  develop  a  grouping  algorithm  for 
near-decomposable  sys  terns . 

We  first  study  the  definition  (3.29).  In  each  area  i,  consisting 

of  n^  machines,  there  are  only  n^-1  independent  functions  z_(t).  We  can 

form  a  basis  of  these  functions  for  the  areas  by  selecting  arbitrarily  one 

machine  which  we  call  reference  machine,  and  forming  the  differences  of  all 

the  other  machines  in  the  same  area  with  respect  to  that  machine.  Doing  so 

for  all  r  areas  we  form  n-r  angle  difference  functions  z^(t).  If  we  denote 

by  x^  the  vector  of  r  machines  selected  as  reference  machines  in  the  areas, 

2 

and  by  x  the  n-r  vector  of  the  remaining  machines,  which  we  will  call 

follower  machines,  then  the  z..  variables  making  vector  z  can  be  written  in 

ij 

terms  of  x  variables  as 


2  T  1 
z  =  x  -Lx  . 


(3.30) 


By  comparing  (3.30)  with  (3.29)  it  can  be  seen  that  the  matrix  L  has  in 

8 

each  row  only  one  nonzero  element  which  is  equal  to  1 .  For  the  row  j , 
corresponding  to  the  j-th  follower  machine,  if  the  element  1  is  in  column  i. 


that  means  that  the  j-th  follower  machine  and  the  i-th  reference  machine  are 

1  2 

in  the  same  area.  This  means  that  given  x  ,  x  ,  and  L  the  coherent  areas 

§ 

are  uniquely  defined.  Therefore  we  call  the  matrix  L  a  grouping  matrix. 
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since  it  groups  the  follower  machines  with  the  reference  machines.  However, 

when  the  areas  are  specified,  there  is  not  a  unique  choice  of  variables  for 
12  12 

x  and  x  ;  to  each  choice  of  x  and  x  there  corresponds  a  different  L^.  As 
an  illustration  consider  a  three  area  five  machine  system.  Given 


1  _  , 

x  =  (x1,x2,x3) 

2  _  ,  >T 
x  =  (x3,x5) 


~1  0  0 

.0  10 


(3.31) 


the  three  areas,  which  are  composed  of  machines  1  and  3,  machines  2  and  5 


ft 


and  machine  4,  are  uniquely  defined.  For  the  same  areas  a  different  choice  of 


1  r  \T 

x  =  (x4,x3,x2) 

x2  =  (x3,x5) ' 


(3.32) 


will  result  in  a  different  L 


0  10 


0  0  1 


(3.33) 


Note  the  zero  column  in  of  (3.31)  or  (3.33)  indicates  the  presence  of  a 
single  machine  area.  Equation  (3.30)  can  be  interpreted  as  a  special  case 
of  the  transformation  (2.11) 


I  0  x  [  x 

=  T 

T  T  2  L  2  ■ 

-L  I  x  x 


(3.34) 


The  substitution  of  (3.34)  into  (3.11)  results  in 
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h 

z  [R(L) 


A12  * 


(3.35) 


where  B^,  and  the  Riccati  equation  R(L)  are  given  by  (2.9),  (2.10),  and 
(2.13).  To  meet  the  definition  (3.29)  system  (3.35)  must  have  the  following 

c 

properties.  First  z  has  to  contain  only  modes  corresponding  to  a^,  for  any 
initial  condition  x(0) .  That  implies  that  R(L) =  0  with  L  being  the  solution 
corresponding  to  the  spectrum  o^.  Next  from  (3.34)  and  (3.29)  it  follows 
that  L  =  L^.  These  observations  allow  an  alternative  formulation  of  coherency 
in  terms  of  the  generalized  Riccati  equation. 

Theorem  3.7;  In  an  ri-machine  system  let  x^  be  the  angles  of  r  noncoherent 
2 

machines  and  x  be  the  angles  of  the  other  n-r  machines.  This  system  is 
r-decomposable  if  and  only  if  the  solution  L  of  R(L)  =  0  corresponding  to  a 
given  is  a  grouping  matrix  L  =  Lg. 

Proof :  According  to  Theorem  2.3,  z  will  have  only  the  modes  corresponding 


to  a  if  R(L)  =  0  such  that  I'm 
r 


is  the  invariant  subspace  of  A  corresponding 


to  o  .  For  the  definition  (3.29)  to  be  satisfied  this  L  must  be  equal  to  an 

L  .  Conversely,  if  a  solution  of  R(L) = 0  exists  and  is  L  ,  it  satisfies  the 
g  g 

definition  (3.29). 

Suppose  now  that  the  system  is  r-decomposable  with  respect  to  a  given  o  ,  but 

we  do  not  know  its  areas.  How  can  Theorem  3.3help  us  find  them?  First,  we 
1  2 

make  a  choice  of  x  and  x  ,  which  in  turn  defines  the  corresponding  equation 
R(L)  =  0.  If  our  x1  does  not  contain  coherent  machines  this  equation  will 
have  the  solution  L  which  is  the  grouping  matrix  needed  to  define  areas.  If 
our  x^  contains  coherent  machines,  the  solution  L  will  not  exist.  The 
negative  outcome  would  mean  that  a  new  choice  of  x^  would  have  to  be  made 


and  a  new  equation  R(L)  = 0  is  solved. 
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The  situation  is  considerably  more  involved  if  this  direct  search 
is  applied  to  a  near-decomposable  system  due  to  the  fact  that  L  is  no  longer 
a  grouping  matrix.  A  further  difficulty  is  that  for  any  choice  of  machine 
angles  x\  the  solution  L  usually  exists.  In  principle,  one  can  make  all 
reasonable  choices  of  x^  and  form  the  set  D  defined  next. 


Definition  3.8:  Given  a  and  a  basis  V  of  the  corresponding  invariant 

r  r.'.  i 


subspace  of  A.  Consider  V =  PV K 


[Yil 

LV2  J  ’ 


where  P  is  any  permutation  matrix 


Define  D  =  {L  =  V2V^}  for  all  P  such  that  exists. 

If  the  system  is  indeed  near-decomposable,  the  set  D  will  contain  at 

least  one  element  sufficiently  close  to  a  grouping  matrix  L  .  One  possible 

S 

measure  of  the  approximation  between  L  and  L  is  a  norm  of  A = L-L  .  The  norm 


is  defined  as  the  maximum  sum  of  absolute  values  of  row  elements,  that  is 


|L||  =  max  E  U.  .  ,  L  =  (£  )  . 

i  j=l 


(3.36) 


Thus  if  there  exists  a  solution  in  D  which  is  also  a  grouping  matrix,  that 

is,  a  decomposable  system,  then  the  norm  (L-L  )  is  zero.  In  a  near- 

g 

decomposable  system,  one  can  search  over  all  the  solutions  in  D  and  find  one 
which  minimizes  norm  (L-L^) .  The  systematic  method  for  area  identification 
in  the  next  section  avoids  this  exhaustive  search.  However,  it  is  interesting 
to  show  an  example  of  the  search  procedure.  Let  us  decompose  into  two  areas 
the  simple  three-machine  system  (Figure  3.2)  whose  model  is 


-14.3  5.4: 

14.3  -49.4 


35.1  x. 


58.4  81.4  -140. 


(3.36) 


Let  the  spectrum  a ^  consist  of  the  two  eigenvalues  of  smallest  magnitude. 
There  are  only  three  possible  choices  of  reference  machines  for  this  example. 


t 


0.723  pu  0.2703pu  D(pu)  =  9.6 

F? -6963 


Line  # 

From 

To 

R(pu) 

X(pu) 

B/2(pu) 

1 

1 

4 

0. 

0.0567 

0. 

2 

4 

5 

0.017 

0.092 

0.079 

3 

5 

6 

0.039 

0.170 

0.179 

4 

3 

6 

0. 

1.0586 

0. 

5 

6 

7 

0.0119 

0.1008 

0.1045 

6 

7 

8 

0.0085 

0.072 

0.0745 

“* 

8 

2 

0. 

0.0625 

0. 

8 

8 

9 

0.032 

0.161 

0.153 

9 

9 

4 

0.01 

0.085 

0.088 

Figure 

3.2.  An 

example  of  the 

3  machine 

system. 

x  *  (x^x^j)  ,  x  =  (x2#x3)  ,  and  x  =  (x^x^,)  and  two  possible  choices  of  Lg, 
[0  1]  and  [1  0].  As  our  first  choice  of  reference  machines  consider 
x"*"  =  (x^.x^)'.  Then  the  solution  of  the  Riccati  equation,  which  is  for  the 
given  spectrum  the  dichotomic  solution  is 


L,  =  1-0.47  1.47]. 

a 


(3.37) 


Of  the  two  possible  L  matrices,  the  one  with  the  best  approximation  of 

g 

L .  is  L  -  [0  1] ,  with 

d  g 


IL.-L  It  =  0.94. 
d  g 


(3.38) 


1  T 

The  second  choice  of  x  =  (x^x^)  will  result  in  a  dichotomic  solution  of  the 
corresponding  equation  R(L) =  0 


L.  *  [-2.13  3.13]. 

d 


(3.39) 


Of  the  two  possible  grouping  matrices  the  one  which  is  a  better  approximation 

of  L,  is  L  =  [0  1],  with 

d  k 


IL.-L  II  =  4.26. 
d  g 


(3.40) 


1  T 

The  third  possible  choice  of  x  =  (x^,,x2)  will  result  in  the  dichotomic  solu¬ 
tion  of  the  corresponding  Riccati  equation 


L  ,  =  [0.32  0.68] 

d 


(3.41) 


which  is  best  approximated  with  L  =  [0  1].  This  gives 

g 


IL.-L  II  =  0.64. 
d  g 


(3.42) 


The  third  choice  of  reference  machines  results  in  the  solution  L  which  can 
be  better  approximated  by  an  L  than  for  any  other  choice.  The  areas  defined 
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by  the  resulting  are  machine  1  in  one  area  and  machines  2  and  3  in  the 
other  area.  There  are  several  other  observations  that  can  be  drawn  from 
this  example.  First,  by  comparing  (3.41)  and  (3.37)  it  can  be  seen  that 
there  is  more  than  one  choice  of  reference  machines  which  result  in  the  same 
area  definition,  although  the  grouping  is  most  clearly  displayed  in  the  third 
case.  The  second  choice  of  reference  machines  results  in  with  the  largest 
norm.  This  is  the  case  when  the  reference  machines  are  taken  from  the  same 
area.  If  this  system  was  r-decomposable ,  we  would  not  have  had  the  solution 
at  all. 

In  remark  2.8  it  was  pointed  out  that  the  spectrum  separation  and 

the  modal  method  order  reduction  are  equivalent  problems.  In  Theorem  3.7  it 

was  shown  that  r-decomposability  can  be  expressed  as  the  spectrum  separation 

problem.  Hence  for  r-decomposable  systems  all  the  three  methods  for  reduced 

order  modeling:  modal  method,  spectrum  separation,  and  coherency  are  directly 

related.  After  the  aggregated  model  based  on  coherent  areas  is  presented  in 

Chapter  4,  it  will  be  shown  that  the  modal  method  and  this  aggregated  model 

based  on  coherency  have  the  same  eigenvalues.  Therefore  one  of  the  tests  for 

validity  of  areas  can  be  the  comparison  of  eigenvalues  of  and  (3.35) 

with  those  of  A  where  and  B^  are  computed  according  to  (2.9)  and  (2.10)  by 

using  L  for  L.  For  our  three  machine  system  o(B,)  =  {0,  -28.6},  o(B0)  =  {-175}, 
8  1  4 

and  ct(A)  =  {0,  -37,  -166},  which  is  another  indication  that  the  areas  are 

near-coherent . 

The  above  direct  search  in  the  three  machine  example  is  presented 
only  as  a  motivation  for  the  systematic  algorithm  presented  in  the  next 


section . 


3.5.  Grouping  Algorithm 


From  the  three  machine  example  it  is  apparent  that  finding  the  areas 

consists  of  two  interdependent  tasks:  first  choosing  the  reference  machines, 

and  second  associating  the  other  machines  to  the  reference  machines.  The 

approach  used  in  the  three  machine  example  is  to  exhaust  all  possible  choices 

of  x^  and  L  ,  that  is  for  each  L,  a  particular  L  was  found  to  minimize  norm 
§  8 

(L-L  ) .  The  best  choice  of  x  is  the  one  corresponding  to  the  smallest  of 
S 

these  minima.  When  the  order  of  the  system  is  large,  this  exhaustive  search 

would  be  computationally  prohibitive.  Due  to  the  properties  of  the  set  D 

established  next  in  Lemmas  3.5  and  3.6,  the  exhaustive  search  can  be  avoided. 

The  algorithm  presented  in  this  section  computes  only  one  element  of  the  set 

D,  which  does  not  necessarily  minimize  norm(L-L  ) ,  but  still  unambiguously 

S 

determines  the  areas . 

Lemma  3.9:  Given  such  that  Oe  o^.  Then  every  element  LeD  defined  by 
Definition  3.4  has  the  property 


(a)  Z  Z,.  1,  i = 1 , 2 , . . . ,n-5 , 

3=1  1J 


(3.43) 


that  is  the  sum  of  row  elements  equals  1,  for  each  row; 


(b)  IILII  >  1. 


(3.44) 


Proof :  From  the  proof  of  Theorem  2.3,  if  u  is  an  eigenvector  of  B^,  then 


(3.45) 


is  an  eigenvector  of  A.  In  particular,  if  then  from  (PI) 
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Thus  (3.43)  is  obtained  by  writing  L,u  in  scalar  form.  The  (b)  follows 

a  o 

from  (3.43)  and  the  definition  of  the  norm  (3.36). 

Notice  that  every  has  the  norm  equal  to  1.  Hence  the  search  for 

a  grouping  matrix  due  to  (3.44)  has  to  be  based  on  the  elements  of  D  with  the 

small  norm.  Based  on  (3.44)  we  can  define  a  grouping  error  as  follows. 

Definition  3.10:  Let  L  be  an  L£  D  with  the  smallest  norm.  Then  a  grouping 
-  m 

error  for  a  given  L  is 

g 

E  =  II L  -L  II  .  (3.47) 

m  g 

From  this  definition  it  follows  that  0<E<°°,  with  E=0  if  and  only  if  the 
system  is  r-decomposable . 

Using  the  connection  between  coherency  and  spectrum  separation 
analyzed  in  the  previous  chapter,  we  can  now  analyze  why  some  solutions  to  L 
have  large  norm.  L  is  computed  using  (2.20)  for  a  given  basis  matrix  V. 

From  the  fact  that  ^  j  is  also  a  valid  basis  of  the  invariant  subspace 
corresponding  to  o^,  it  follows  that  any  other  basis  will  have  identically 
the  same  row  vectors  in  V  corresponding  to  the  machines  in  the  same  area. 
Geometrically,  there  will  be  r  distinct  groups  of  identical  row  vectors  of  V. 
In  the  case  of  near-coherency  we  will  have,  instead,  r  narrow  nonoverlapping 
cones  containing  all  the  row  vectors  of  V.  Now  it  is  clear  that  if  reference 
machines,  which  give  components  of  the  vector  x^,  contain  two  machines  from 
the  same  area,  then  the  resulting  will  have  two  identical  row  vectors  (in 
the  case  of  r-decomposable  systems,  or  near-identical  in  t  case  of  r-near- 
decomposable  systems)  and  hence  the  solution  will  not  exist,  or  will  have  a 


very  large  norm.  That  was  the  case  with  the  second  selection  of  reference 
machines  in  the  example  of  the  previous  section.  Since  at  the  time  of 
computing  a  basis  V  it  is  generally  not  known  whether  the  selection  of 
reference  machines  is  good  or  not,  the  question  is  whether  the  effort  spent 
on  computing  V  (which  is  considerable)  is  lost.  The  following  result 
answers  the  question. 

Lemma  3.11:  Given  a  subspectrum  of  A  and  a  basis  of  the  corresponding 

invariant  subspace  V.  Then  each  element  L  of  D  is  a  solution  of  the  Riccati 

equation  corresponding  to  for  the  system  with  state  vector  x =  Px . 

Proof :  If  V  spans  an  invariant  subspace  of  A  then  PV  spans  the  invariant 

T 

subspace  corresponding  to  the  same  subspectrum  of  PAP  .  The  result  then 
follows  from  Theorem  2.3.  This  result  shows  that  to  find  an  L  in  D  with  a 
small  norm  we  need  to  compute  a  basis  V  for  a  given  state  ordering,  and  then 
perform  a  suitable  permutation  of  rows.  In  finding  such  a  permutation  we 
benefit  from  the  geometrical  interpretation  of  coherency  requiring  existence 
of  narrow  cones  of  row  vectors  in  V.  Hence,  to  ensure  invertibility  of 
we  want  to  select  r  vectors  from  the  r  different  cones.  In  this  way  we  ensure 
invertibility  and  well  conditioning  of  V^.  To  find  this  set  of  r  rows,  we  use 
Gaussian  elimination  with  complete  pivoting.  During  the  elimination,  the  rows 
and  columns  of  V  are  permuted  such  that  the  (1,1)  entry  of  the  resulting  V  is 
the  largest  entry  in  magnitude.  Note  that  permuting  the  rows  of  V  is  equi¬ 
valent  to  changing  the  ordering  of  the  machines.  This  (1,1)  entry  of  V  is 
used  as  the  pivot  for  performing  the  first  step  of  the  Gaussian  elimination. 
Then  the  largest  entry  is  chosen  from  the  remaining  (n-l)x(r-l)  submatrix  of 
the  reduced  V  and  is  used  as  the  pivot  for  the  next  elimination  step.  The 
elimination  terminates  in  r  steps  and  the  machines  corresponding  to  the  first 
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f 

B 


r  rows  of  the  final  reduced  V  matrix  are  designated  as  the  reference  machines. 
In  this  Gaussian  elimination  process,  rows  having  small  entries  will  not  be 
used  as  the  pivoting  row  because  these  small  entries  are  the  result  of  eli¬ 
mination  with  almost  identical  rows  already  used  as  pivoting  rows.  Thus, 
this  algorithm  does  not  put  two  near-coherent  machines  together  as  reference 
machines . 

For  the  set  of  reference  machines  found  by  the  algorithm  the 
corresponding  L  is  readily  computed  from 

T  T  T 

V^L1  =  Vj  (3.48) 

using  the  LU  decomposition  of  obtained  from  the  Gaussian  elimination.  The 

next  step  is  to  find  an  L  approximating  L,  that  is  to  find  the  machines 

§ 

belonging  to  each  area. 

From  (3.47)  we  see  that  if  row  vectors  of  are  understood  as  a 
basis  for  an  r-dimensional  space,  then  elements  of  L  are  coordinates  of  the 
row  vectors  of  in  this  basis.  Specifically,  the  elements  of  a  j-th  row 
of  L,  are  coordinates  of  the  j-th  vector  in  on 

Vf ,vz , . . . ,vr •  If  we  now  recall  the  geometrical  image  of  coherent  areas  as 
consisting  of  narrow  nonoverlapping  cones  of  row  vectors  in  V,  and  that  the 
Gaussian  elimination  selects  r  vectors  from  r  different  cones,  it  becomes 
clear  that  in  each  row  j  of  L,  corresponding  to  a  machine  which  belongs  to 
some  area  k,  there  will  be  only  one  element  close  to  1,  which  is  the  pro¬ 
jection  of  the  j-th  vector  of  on  the  basis  vector  v^  along  the  space 
spanned  by  vectors  v^,  i=l,2,...,r,  i/k;  all  the  other  coordinates  of  this 
vector,  i.e.,  all  the  other  elements  in  the  same  row  of  L  will  be  small. 


72 


Therefore  to  find  an  L  from  a  computed  L  we  proceed  as  follows.  In  each 

s 

row  of  L  we  find  the  largest  positive  element,  approximate  it  by  one  and 

approximate  all  other  elements  in  the  same  row  by  zero.  In  view  of  the  above 

analysis  this  is  a  suboptimal  procedure  for  minimizing  the  norm  of  L-L  ,  i.e., 

8 

the  grouping  error  of  Definition  3.10.  Although  we  cannot  always  claim  that 

II LB  <  9 L T II  ,  L'€D  (therefore  the  suboptimality  of  the  algorithm),  it  is  clear 

from  the  geometrical  explanation  that  for  near  r-decomposable  systems  the 

algorithm  gives  a  unique  definition  of  areas.  In  other  words  if  the 

"coherency  cones"  are  sufficiently  narrow,  no  matter  which  reference 

vectors  are  selected  as  long  as  there  are  r  vectors  from  r  different  cones, 

the  belonging  of  the  other  vectors  to  these  cones  is  uniquely  defined  by  the 

largest  element  in  each  row  of  L.  In  this  case  II L-L  II  =  0(e)  ,  where  e  is  the 

S 

largest  angle  among  two  angles  in  the  same  cone.  As  the  angles  of  the  cones 
start  increasing,  the  criterion  for  machine  grouping  based  on  the  largest 
element  in  each  row  of  L  becomes  less  discriminative.  For  each  follower 
machine  we  can  define  a  discrimination  factor  as  follows. 

Definition  3.12:  Given  an  LGD  with  the  smallest  norm,  let  i  be  the  column 

index  of  the  largest  positive  element  and  k  be  the  column  index  of  the  second 

largest  positive  element  in  the  same  row.  The  discrimination  factor  for 

the  follower  machine  j,  DF.  is 

J 


DF.  =  III  £..-*..11  -  IU.-i.  II 
J  J  i  J  k 


(3.49) 


where  l.  is  the  i-th  row  of  L  and  e.  is  the  i-th  row  of  the  (r*r)  identity 
J  i 


matrix. 


The  largest  value  of  DF  is  2,  which  is  in  the  case  of  misplacing 


ideally  coherent  machines,  and  the  smallest  value  is  zero,  which  is  in  the 


case  when  v.ev-  has  the  same  coordinates  on  both  v.  and  v,  ,  v.,v,  €  v,  . 

j  2  1  k  1  k  1 

Machines  with  the  small  discrimination  factor  should  be  given  special 
consideration.  First,  they  can  be  grouped  based  on  approximating  the  second 
largest  element  in  the  row  by  one  if  that  choice  is  more  convenient  for  some 
reason  (for  example  to  achieve  an  area  grouping  consistent  with  the  administra¬ 
tive  boundaries  of  the  areas).  Second,  in  the  case  of  network  changes  these 
are  the  machines  most  likely  to  change  their  area  belonging. 

We  now  summarize  the  grouping  algorithm  as  follows: 

Step  1:  Decide  on  the  number  of  areas  and  spectrum  a  . 

Step  2:  Compute  a  basis  matrix  V  for  a  given  ordering  of  the  x  variables 

and  a  . 
r 

Step  3:  Apply  Gaussian  elimination  with  complete  pivoting  to  V  and  obtain 
the  set  of  reference  machines. 

Step  4:  Compute  L  for  the  set  of  reference  machines  chosen  in  step  3, 

Step  5:  Construct  the  matrix  L  ana  find  the  machines  in  each  area. 

g 

The  main  computational  load  is  in  step  2.  However,  only  a  partial  eigen- 
subspace  V  of  A  is  required  and  since  A  is  similar  to  a  symmetric  matrix, 
eigenvalue-eigenvector  computation  is  well  conditioned  [46].  The  numerical 
aspects  of  the  basis  computation  will  be  discussed  in  more  detail  in 
Section  4.8. 

We  illustrate  this  area  selection  procedure  on  a  48  machine  system 
from  [18].  The  data  are  given  in  the  reference.  From  the  linearized  system 
state  matrix  A,  we  extract  the  n*n  matrix  A  of  (3.11)  by  eliminating  rows  of 
A  corresponding  to  <5  and  columns  corresponding  to  w.  This  matrix  is  given 
in  Appendix  B.  In  the  first  step  we  specify  that  we  want  nine  areas  with 


respect  to  the  nine  slowest  modes.  From  this  point  on  the  algorithm  proceeds 
automatically  giving  the  following  results.  In  step  2  a  basis  for  the  9- 
dimensional  slow  subspace  is  computed.  In  step  3  the  Gaussian  elimination 
is  performed  and  the  set  of  reference  machines  is  found  to  be  5,  39,  44,  34, 
48,  41,  17,  29,  36.  In  Step  4  the  solution  L  is  found  (for  this  selection  of 
the  spectrum  it  will  be  a  dichotomic  solution)  and  it  is  given  in  Table  3.1. 
The  largest  element  in  each  row,  those  elements  are  underlined  in  Table  3.1, 
is  used  to  identify  the  belonging  of  the  corresponding  machine  to  an  area, 
i.e.,  this  element  is  approximated  by  1  and  all  other  elements  in  the  row  by 
zero.  As  a  result  the  following  grouping  of  machines  into  areas  is  obtained: 
Area  1:  1,  2,  3,  4,  5,  6,  7,  8,  9 

Area  2:  39,  42 

Area  3:  43,  44,  45,  46 

Area  4:  34,  35 

Area  5:  48 

Area  6:  32,  37,  38,  40,  41 

Area  7:  13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  31,  47 

Area  8:  10,  27,  28,  29,  30 

Area  9:  11,  12,  33,  36. 

The  grouping  error  (3.47)  normalized  with  the  number  of  follower  machines  is 
0.74.  The  average  value  of  the  elements  of  L  approximated  by  one  is  .63  and 
the  average  value  of  elements  of  L  approximated  by  zero  is  .05.  In  average, 
the  discrimination  factor  is  0.85  (compared  with  the  maximal  value  of  2). 
However,  for  some  machines  like  10,  14,  24,  25,  31,  32,  33,  37,  38,  42,  and  47 
the  discrimination  factor  is  much  smaller  than  the  average  value.  For  example 


75 


Table  3.1.  Ld  for  the  48  machine  system. 


£ 


D 


0 


5 

39 

44 

34 

48 

41 

17 

29 

36 

1. 

0. 56 

-0.00 

-0.04 

-0. 01 

-0. 00 

-0. 00 

0.  05 

0. 42 

0.02 

2. 

oTSo 

-0. 00 

-0.04 

-0. 01 

-0. 00 

-0.00 

0. 05 

0.  38 

0.02 

3. 

0.83 

-0.00 

-0.01 

-0. 01 

0.  00 

-0. 00 

0. 03 

0.  15 

0.01 

4. 

0.83 

-0. 00 

-0. 01 

-0.00 

0. 00 

-0. 00 

0.  03 

0.  14 

0. 01 

6. 

0. 84 

-0. 00 

-0.02 

-0. 01 

-0. 00 

-0.00 

0.03 

0.  15 

0. 01 

7. 

0.  83 

-0. 00 

-0.02 

-0. 01 

-0.00 

-0.00 

0. 03 

0.  16 

0. 01 

8. 

0. 85 

-0. 00 

-0. 01 

-0. 01 

0. 00 

-0. 00 

0.05 

0.  11 

0.  02 

9. 

0.58 

-0. 00 

-0. 01 

-0. 01 

0.00 

-0. 01 

0.  10 

0.  30 

0. 04 

10. 

6776 

-0.00 

0.03 

-0.03 

0. 01 

-0.01 

0.28 

0.43 

0.  13 

11  . 

0.  07 

-0. 01 

-0.03 

-0.  15 

0.00 

-0.04 

0. 21 

0. 28 

0. 66 

12. 

0.  08 

-0. 01 

-0. 02 

-0.  12 

0.  00 

-0. 03 

0. 21 

0. 27 

0. 62 

13. 

0.  10 

-0. 00 

-0. 01 

-0.03 

0.00 

-0. 02 

0.  51 

0. 32 

0.  12 

14. 

0.  11 

-0. 00 

0.02 

-0. 02 

0. 01 

-0. 01 

0. 46 

0.  32 

0.  12 

15. 

0.04 

-0.00 

-0. 01 

-0.02 

-0.00 

-0.02 

0.78 

0.  17 

0. 06 

16. 

0.  02 

0. 01 

-0.00 

0. 03 

0.00 

0. 01 

0.77 

0.  10 

0.  06 

18. 

0. 02 

0.  00 

0.  00 

0.02 

0. 01 

0.  01 

0.73 

0.  09 

0. 06 

19. 

0.03 

0. 01 

0.  03 

0.04 

0.  02 

0.02 

0. 67 

0.  10 

0.  08 

20. 

0.  03 

0. 01 

0.  02 

0.  03 

0.02 

0.  01 

0.72 

0.  10 

0. 06 

21  . 

0.  02 

0.00 

0.  11 

0. 02 

0. 08 

0.  00 

0. 63 

0. 09 

0. 05 

22. 

0.  03 

0. 00 

0. 09 

0.  02 

0.  07 

0. 01 

7T75T 

0.  10 

0. 06 

23. 

0.06 

0. 00 

0. 05 

-0. 00 

0.  01 

-0. 01 

0762 

0.  19 

0.  07 

24. 

0.  09 

0.  00 

0.  21 

-0.01 

0.  04 

-0.00 

0.  37 

0.24 

0.06 

25. 

0.  11 

-0. 00 

0.  17 

-0. 02 

0.03 

-0. 01 

6736 

0.  30 

0. 07 

26. 

0.  10 

-0. 00 

0. 09 

-0.03 

0. 02 

-0.02 

0.  42 

0.  33 

0. 09 

27. 

0. 07 

-0.00 

0.  12 

-0.00 

0. 01 

-0.00 

0.05 

0.73 

0.  02 

28. 

0.02 

-0. 00 

0. 03 

-0.00 

0.00 

-0.00 

0. 01 

0.93 

0.00 

30. 

0.  04 

-0. 00 

0.06 

-0.00 

0.00 

-0.00 

0.  02 

0.  88 

0. 01 

31  • 

0.  02 

0.05 

0.  00 

0.  13 

0. 01 

0.  10 

0.  48 

0. 06 

0.  15 

32. 

0.  00 

0. 24 

-0.00 

0.09 

0.  04 

0.  41 

0.  16 

-0.01 

0. 07 

33- 

0.  00 

0.04 

-0.00 

0. 23 

0. 00 

0.  10 

0. 25 

0. 02 

0.  31 

35. 

0. 00 

0. 01 

-0.00 

0. 87 

-0.00 

0.02 

0.05 

0. 00 

0.05 

37. 

0.  00 

0.  32 

0.01 

0. 05 

0.09 

0. 40 

0. 09 

-0. 01 

0. 04 

38. 

0.  00 

0. 38 

0.00 

0. 04 

0.06 

0. 41 

0. 03 

-0. 01 

0. 03 

40. 

0.  00 

0. 28 

-0.00 

0. 05 

0. 00 

0.  56 

0. 09 

-0. 01 

0.  03 

42. 

0.  00 

0. 47 

0.00 

0.03 

0. 03 

0.  39 

0. 06 

-0. 01 

0.  02 

43- 

0. 01 

0. 00 

0.73 

-0. 00 

0.  1 1 

0. 00 

0.  02 

0.  13 

0.  00 

45. 

0.  02 

-0. 00 

676o 

-0.00 

0.  12 

-0.00 

0.  18 

0. 07 

0. 02 

46. 

0.  00 

0. 00 

0. 89 

0. 00 

0. 07 

0. 01 

0. 01 

0.  01 

0.  00 

47. 

0.02 

0. 00 

0. 26 

0. 01 

0. 21 

0. 00 

0. 08 

0.  04 
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k 
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for  machine  33,  DF^  =  0.06.  All  the  machines  with  a  small  discrimination 
factor  will  be  called  critical  machines.  Note  from  Figure  3.3  that  all  the 
critical  machines  are  border  machines  of  their  respective  areas,  but  not  all 
the  border  machines  are  critical  machines.  Despite  the  presence  of  the 
critical  machines,  the  area  grouping  is  good,  as  will  be  shown  in  Section  4.7 
by  some  characteristic  machine  responses. 


3.6.  Area  Variables  and  Intermachine  Variables 

In  the  previous  section  it  was  shown  that  the  coherency  is  equiva¬ 
lent  to  the  requirement  that  a  grouping  matrix  achieves  the  spectrum  separation. 
Then,  the  state  variables  associated  with  the  subsystem  with  state  matrix 
are  z  variables  (3.29),  i.e.,  intermachine  variables  within  coherent  areas. 

X 

In  this  formulation  the  subsystem  (x  ,B^)  is  still  coupled  to  the  subsystem 
(zjB^).  For  the  separate  analysis  of  these  two  subsystems  it  is  necessary  to 
decouple  them,  and  one  way  to  do  it  is  by  using  the  transformation  (2.52) 
introduced  in  Section  2.4, 


r 


y 

z 


I 

0 


(3.50) 


Applying  this  transformation  to  (3.35)  with  L  satisfying  R(L)  = 0  results  in 


y 

z 


(3.51) 


As  shown  in  Section  2.4,  equation  P(H)  = 0  has  always  the  solution,  whenever 

Q 

c  <"i  cr  =0.  However,  in  this  particular  case  even  more  can  be  concluded, 
r  r 
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LgJ5na_3.13:  Consider  the  matrix  A=-M~1K  of  (3.11),  where  K  is  symmetric  and 
M  is  the  diagonal  matrix  of  machine  moments  of  intertia,  whose  (rxr)  and 
(n-r)x(n-r)  diagonal  blocks  are  and  M2,  respectively.  Then  each  solution 
L  of  R(L) = 0  corresponding  to  some  o^,  and  the  corresponding  solution  H  of 
P(H)  =  0  are  related  by 


H  =  (M^  +  L,TM2L)-1LTM2  . 


(3.52) 


Proof:  Let  us  first  rewrite  R(L) *  0  and  P(M)  =  0  as 


c-l  n*  l-o 


[I  H]  B  =  0. 

I 


(3.53) 


(3.54) 


Substituting  B  =  T^AT^  into  (3.54)  yields 


[I-HL  H]  A  =  [I-HL  H](-M~1K)M  1M 


-1  -I  F  "MtH 

-  [(I-HL)M  1  HM  L]  AT  1  =  0. 

|_  M  (  T— LH) 


(3.55) 


Pre-  and  post-multiplying  (3.53)  by  [M2(I-LH)]T  and  [(I-HLJM"1]1  and 
comparing  to  the  transpose  of  (3.55),  we  obtain 

H1^  =  (I-LH)TM2L 


(3.56) 


which  simplifies  to  (3.52). 


The  same  relationship  can  be  obtained  by  modal  methods  [16J.  Under 
the  conditions  of  this  lemma  the  complete  transformation  from  x  to  (y.z) 


« 
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variables  is 


M  1M1  M_1LTM„ 

a  1  a  2 

-L  I 


(3.57) 


where 

M  =  M.  +  LTM-L.  (3.58) 

a  1  2 

If  the  system  is  r-decomposable  then  H  has  the  same  "zero-nonzero"  pattern 
T 

as  L  and  we  denote  it  by  H  .  With  L  and  H  used  in  (3.57),  an  interesting 
g  g  g  g 

interpretation  of  the  physical  meaning  of  y  and  z  variables  follows.  First, 

notice  from  (3.58)  that  in  this  case  M  is  the  matrix  of  area  inertias,  i.e. 

a 

the  i-th  entry  contains  the  sum  of  moments  of  inertia  of  all  the  machines  in 
the  i-th  area.  Then  from  the  first  row  of  (3.57)  we  get  immediately  that 
y  variables  are  familiar  "area  center  of  inertia"  [41],  given  by 


y.  =  l  M.x./M  .,  for  all  i  in  area  i. 
i  j  J  J  ax’  J 


(3.59) 


In  summary,  in  the  case  of  r-decomposable  systems  by  dividing  the  system  into 
coherent  areas,  it  is  possible  to  construct  directly  two  low  order  subsystems, 
one  describing  the  dynamics  of  intermachine  oscillations  within  coherent  areas 
and  the  other  describing  the  interarea  dynamics.  For  the  r-near-decomposable 
system  we  still  apply  the  same  transformation  (3.59)  and  (3.30).  We  will  get 
the  two  subsystems  which  are  weakly  coupled,  rather  than  completely  decoupled. 
The  same  is  true  for  models  with  damping  and  nonlinear  models.  When  the 
spectrum  defining  coherency  is  given  as  r  slowest  eigenvalues  of  A,  the  pro¬ 
cedure  will  result  in  y  variables  being  the  slow  states  and  z  variables  being 
the  fast  states  of  the  system.  For  the  r-near-uecomposable  systems,  their 
basic  slow,  i.e.,  fast  character  will  be  retained,  which  means  that  by  means 


of  coherency  the  mixed  state  model  (3.11)  is  transformed  into  singularly 
perturbed  form.  In  [25]  it  was  shown  that  application  of  singular 
perturbation  on  the  transformed  model  gives  excellent  results  in  terms  of 
eigenvalue  approximation  of  the  reduced  subsystems  compared  with  the  full 
system  eigenvalues.  More  properties  of  the  slow  coherency  will  be  given  in 
the  next  chapter. 


4.  SLOW  COHERENCY,  WEAK  COUPLING,  AND  NONLINEAR  EQUIVALENTS 

4.1.  Introduction 

In  this  chapter  we  show  that  only  if  coherency  is  defined  with 
respect  to  the  slowest  modes  of  the  linearized  system  (3.11),  the  resulting 
areas  are  weakly  coupled,  i.e.,  weaker  than  for  any  other  partial  coherency. 

The  area  boundaries  are  in  this  case  relatively  insensitive  to  small  loads 
and  paramete:  changes.  These  two  factors  taken  together  indicate  certain 
robustness  of  area  boundaries,  making  them  useful  for  nonlinear  system 
analysis.  Possibilities  for  using  slow  coherency  for  analysis  of  nonlinear 
models  are  discussed  and  conclusions  verified  on  a  48  machine  system.  Benefits 
in  numerical  identification  of  coherent  areas  resulting  from  using  slow 
coherency  are  discussed  at  the  end. 

There  is  no  previous  published  work  on  the  relation  between 
coherenc)  and  weak  coupling,  as  well  as  the  analytical  study  of  the  sensiti¬ 
vity  of  area  boundaries  to  parameter  changes.  However,  the  use  of  linear 
coherency  for  aggregation  of  nonlinear  models  is  implicit,  or  alluded  to  in 
[16,19,41,43].  Other  uses  of  slow  coherency  considered  in  this  chapter  [45], 
which  exploit  the  connection  between  time  scales  and  coherency,  are  based  on 
singular  perturbation  theory. 

4.2.  Electromechanical  Equivalent  of  Area  Dynamics 

As  a  preparation  for  exposition  of  the  main  result  of  this  chapter 
concerning  a  relation  between  slow  coherency  and  weak  coupling  between  areas, 
we  first  study  aggregated  models  of  the  linearized  system  (3.11).  This 


analysis  is  extended  to  nonlinear  models  in  Section  4.6.  Once  areas  are 
identified,  their  aggregation  is  intuitively  appealing  to  a  power  system 
expert  [19,41],  Each  area  is  represented  by  one  generator,  those  are  then 
suitably  interconnected  to  preserve  the  total  power  interchange  between 
areas.  The  resulting  interconnection  matrix  is  generally  nonsymmetric  [19]. 
In  this  section  we  proceed  differently.  We  use  analytical  tools  of  Chapter  3 
to  establish  aggregated  models  and  their  properties.  One  of  the  results  is 
that  the  aggregated  model  has  the  same  structure  as  in  previous  approaches, 
but  the  interconnection  matrix  is  symmetric,  so  that  it  can  be  modeled  as  a 
passive  network  without  phase  shifters.  The  result  is  given  in  the  following 
lemma . 

Lemma  4.1;  Let  the  coherency  in  the  system  (3.11)  be  defined  with  respect 

to  a  spectrum  a  containing  zero  eigenvalue,  and  let  L  :  R(L  ) = 0  define 
^  6  § 

corresponding  areas.  Then 

(i)  the  aggregated  model  of  (3.11)  which  reproduces  area  dynamics 
is  defined  by 


y  =  Sx 


(4.1) 


where 


S  =  (I-H  L 
g  g 


+  H  )  . 
g 


(4.2) 


(ii)  The  aggregated  state  matrix  F  =  can  ^actoret^ 


as 


B,  =  M  LK 
1  a  a 


(4.3) 


where  M  is  the  diagonal  matrix  of  area  inertias  (3.50)  and  K  is  a  symmetric 
a  a 


matrix  defined  as 
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K  =  (I  L  )  K 
a  g 


(4.5) 


(iii)  is  a  function  of  connections  between  machines  in 
different  areas  only  and  not  a  function  of  connections  between  machines  in 
the  same  area,  i.e. 


k.  .  = 

ija 


E  E  k  ,  i^j 

pel(i)  q€l(j)  pq 


r 

-E  k. 


(4.5) 


Before  giving  the  proof  let  us  emphasize  the  following.  In  the  standard 

aggregation  approaches  [6,7],  the  main  concern  is  to  obtain  a  low  order  model 

which  reproduces  a  given  spectrum.  In  application  to  the  electromechanical 

model  of  power  systems,  once  L  and  the  corresponding  H  are  known,  more  can 

g 

be  achieved.  Namely,  the  complete  decomposition  of  the  system  into  two  sub¬ 
systems  is  possible  in  such  a  way  that  one  subsystem  reproduces  the  given  sub¬ 
spectrum,  and  the  other  subsystem  reproduces  the  complementary  subspectrum. 
Proof :  (i)  The  form  of  S  follows  from  (3.55)  and  (3.50).  In  comparison  to 

the  general  form  of  aggregation  matrix  given  in  [5],  S  =  P(I  0)X  \  it  can  be 
seen  that  in  this  case  P  =  X^,  i.e.,  it  is  the  r-dimensional  principal  submatrix 
of  the  modal  matrix  X. 

(ii)  Again  with  reference  to  the  general  form  of  the  aggregated 


state  matrix  [4],  F =  B^ = SAS  ,  we  see  from  (2.7)  that  the  r.h.  side  inverse 

+  +  T  -IT 

S  is  S  =  (I  L  )  .  Using  the  particular  form  of  the  solution  H  =M  L  M.  (as 
g  g  a  g  2 

-1  T 

given  in  (2.50)),  we  get  S = (I  +  L^)  M  from  which  the  result  follows 
immediately • 
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(iii)  The  result  (4.5)  follows  by  noting  the  specific  "zero-one" 
structure  of  L  and  the  fact  that  the  diagonal  elements  of  K  are  the  negative 

s 

sums  of  off-diagonal  elements.  The  above  properties  are  best  illustrated  by 
an  example. 

Example  4.1;  For  the  (hypothetical)  system  in  Figure  4.1,  the  matrices  M 
and  K  are  as  given  below.  We  will  find  a  grouping  matrix  for  this  system,  and 
then  analyze  the  corresponding  aggregated  model  B^. 


M  *  diag(5 ,1,4, 


K  = 


2 

-3 

1 

0 

0 


0 

1 

-4 

1 

2 


1,1) 

2 

0 

1 

-3 

0 


1 

0 

2 

0 

-3 


If  r=2  and  o9={0,  X^},  then 


L 

g 


0 

1 

1 


satisfies  Riccati  equation  (2.13)  and  defines  coherent  areas  as  indicated  in 
the  figure  by  the  dotted  lines.  The  matrix  H  is  given  by 


-1  T 

H  =  M  L  M- 
g  a  g  2 


4 

9 

0 


0 

1 

3 


0 

1 

3 


Then  from  (2.9)  is 


4.3.  Structural  Conditions  for  Coherency 

Using  the  fact  that  in  the  case  of  coherency  L  and  satisfy 

O 

Riccati  (2.13)  and  Lyapunov  equations  (2.49),  respectively,  allows  alterna¬ 
tive  formulation  of  the  coherency  directly  in  terms  of  network  and  machine 
parameters.  The  following  lemma  states  the  result. 
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Lemma  4.2:  A  necessary  and  sufficient  condition  for  a  system  (3.11)  to  be 


divided  into  r  coherent  areas  according  to  Definition  3.1  is  that 


E  (y A* cos  ek*  .  vr£Vcos  6U  =  0 

tel(i)  \  Mi 


(4.6) 


or  equivalently 


i=l,2,...,r;  j  =  r+l,...,n 


Mi  2,€I(k)Vi  £  i£  i£ 


— —  l  V  V. b  .  cos  6  . 
Mj  ilel(i)  1  2  £j  £l 


i»l,  —  ,r;  j  =  r+l,...,n. 


(4.7) 


For  simplicity  of  notation  it  is  assumed  k  =  k(j)  to  mean  the  repre¬ 
sentative  machine  of  the  area  containing  j. 

Proof :  Condition  (4.6)  follows  from  (2.13)  directly  by  inserting  L  for  L, 

o 

and  (4.7)  follows  from  (2.49)  by  inserting  L  and  H  for  L  and  H. 

g  g 

Similar  conditions  were  also  obtained  in  [21]  using  a  different 
approach.  These  conditions  are  interesting  in  that  they  emphasize  the  role 
of  network  and  machine  parameters  in  forming  areas.  For  most  networks  and 
loading  conditions  cosine  terms  in  (4.6)  and  (4.7)  can  be  approximated  by  1. 
Then  the  two  conditions  say  that  if 


d  =  z  i=l  2 

dji  £si(1)  M.  >  1-2'-’r 

j  =  r+1, . . . ,n 


(4.8) 


is  defined  as  electromechanical  distance  of  machine  j  from  the  area  i,  then 
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« 


[ 


r 

l 

« 


i 


« 


all  the  macines  in  the  coherent  area  k  have  the  same  distance  to  a  given 
area  i,  where  i,k = 1,2, . . . ,r. 

For  the  48  machine  system  of  Section  3.5,  Table  4.1  gives  the 
average  variation  of  this  distance  for  each  of  the  nine  areas.  The  conditions 
of  the  lemma  guarantee  partial  coherency,  but  they  do  not  take  explicitly  into 
account  the  spectrum  corresponding  to  the  coherency. 


4.4.  Weak  Coupling  Between  Areas 

So  far  we  have  considered  partial  coherency  with  respect  to  any 
prespecified  spectrum  a^,  with  the  only  restriction  that  it  had  to  contain  the 
zero  eigenvalue.  In  this  section,  as  well  as  in  Section  4.5,  we  will  show 
that  the  slow  coherency,  i.e.,  the  case  when  contains  r  slowest  eigenvalues 
of  (3.11),  has  particularly  useful  properties. 

The  basis  for  the  following  analysis  is  the  result  of  Section  4.2, 
which  establishes  a  link  between  coherency  defining  spectrum  and  the  strength 
of  connections  between  areas.  Before  stating  the  new  result,  we  introduce 
several  definitions,  similar  to  the  definition  of  the  electromechanical 
distance  (4.8). 

Definition  4.3:  The  electromechanical  connection  between  areas  i  and  j  is 
defined  as 


D. .  =  M-1  Z  Z  k  , 

ai  p€  I(i)  q€I(j)  Pq 


i.j  = 1,2, . . . ,r. 


(4.9) 


4 


Table  4.1.  Average  violation  of  the  criterion  (4.6)  for  the  48  machine 
system. 

Area  ->-123456789 
4- 

1 

2 

3 

4 

5 

6 

7 

8 


.09 

- 

- 

- 

.02 

- 

- 

.04 

- 

- 

.07 

- 

- 

- 

- 

- 

.13 

- 

.02 

- 

.01 

.06 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

- 

.03 

- 

- 

- 

.05 

.02 

.01 

.01 

■ 

.03 

- 

.01 

.01 

.09 

.02 

.01 

.03 

- 

.01 

.  .  . 

- 

- 

- 

.07 

.10 

.02 

- 

.01 

- 

.01 

.02 

.01 

0 
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The  electromechanical  connection  of  machine  i  to  the  system  is 


D,  =  ED.,. 

1  j=l 


(A. 10) 


The  total  interconnection  of  the  system  is 


DT  =  ED.. 
T  i=1  x 


(4.11) 


Assuming  that  the  cos  terms  in  k  are  approximately  1,  then  all  these 
connection  measures  are  basically  measures  of  electromechanical  distance.  We 
also  define 

Definition  4.4:  The  measure  of  the  maximal  strength  of  the  aggregated  model 


T 

XI  X 

a 

m  *  max  — = - 

8  X  xtm  X 
a 


(4.12) 


which  is  the  biggest  ratio  of  the  potential  over  the  kinetic  energy  in  the 
aggregated  model. 

We  are  now  ready  to  characterize  different  types  of  coherency. 

Lemma  4.5:  Suppose  that  in  the  system  (3.11)  there  are  several  different 

L  :R(L  )=0,  each  corresponding  to  a  different  spectrum  a  .  Suppose  that  one 

8  g  r 

contains  r  slowest  eigenvalues  of  the  system,  i.e.,  it  defines  slow 
coherency.  Then, 

(i)  the  total  interconnection  of  the  system,  D^,  is  the  smallest 
for  the  slow  coherency  case; 


(ii)  the  maximal  strength  of  interconnections,  mg ,  is  the  smallest 


for  the  slow  coherency  case. 


Proof :  The  proof  uses  the  relation  between  eigenvalues  of  a and  area  inter¬ 
connections  established  in  Lemma  4.1,  and  various  definitions  of  eigenvalues. 

r 

(i)  This  property  is  established  using  TrB,  =  Z  A. ,  A.ea  ,  the  form  of 

1  i=1  i  i  r 

diagonal  elements  of  which  are  given  by  (4.5),  and  the  fact  that  this  trace 
is  the  smallest  for  the  slow  coherency  case,  (ii)  This  property  can  be 
derived  using  eigenvalue  definition  via  the  Raleigh-quotient  [27]  and  the 
fact  that  the  largest  of  the  slow  eigenvalues,  which  is  equal  to  m,.,  is  the 
smallest  in  the  case  of  slow  coherency.  The  results  of  the  lemma  indicate 
that  slow  coherence  has  the  property  to  divide  the  system  into  weakly  coupled 
areas  so  that  the  total  sum  of  interconnections  is  minimized.  Since  for  a 
given  system  the  total  sum  of  connections  between  machines  is  equal  to  the 
total  sum  of  interconnections  between  areas  and  the  sum  of  connections  inside 


areas,  it  follows  that  the  slow  coherency  also  maximizes  the  total  sum  of 
internal  connections  in  the  areas.  This  property  will  be  further  clarified 
later  on  in  this  section.  With  this  result  an  interpretation  of  the  slow 
coherency  can  be  the  following:  it  is  a  result  of  groups  of  machiens  being 
strongly  coupled  inside  the  groups  and  weakly  coupled  between  the  groups.  The 
results  concerning  the  slow  coherency  are  global,  that  is,  for  the  whole 
system  rather  than  for  each  individual  subsystem,  although  as  we  will  see 
later  on  the  weak  coupling  property  can  be  stated  as  the  average  property  of 
each  subsystem,  which  is  also  observed  on  the  example  system.  Table  4.2  gives 
interconnection  patterns  for  the  48  machine  svstem  divided  into  9  areas  based 
on  slow  coherency.  Entry  t  in  the  table  gives  relative  electromechanical 
connection  between  areas  i  and  j,  where  scaling  is  done  with  the  total  con¬ 


nection  of  area  i,  i.e., 


t. .  =  D. ,/D. . 
ij  ij  i 


Entries  less  than  1%  in  magnitude  are 


denoted  by  -.  Several  observations  can  be  made  from  the  table. 
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D 


(i)  In  all  cases  >  t^  ,  i.e.,  connections  inside  areas  are 
stronger  than  the  connections  between  areas. 

(ii)  The  cases  of  stronger  interconnections  (for  example  "t  ’’) 

2d 

are  the  cases  where  some  machines  in  areas  2  and  6  are  grouped  based  on  rows 

of  L,  which  were  not  close  to  zero-one  pattern,  see  Table  3.1. 
d 

(iii)  The  state  matrix  A  has  the  same  block-diagonal  structure  as 
Table  4.2,  if  the  rows  are  ordered  so  that  all  the  machines  in  area  1  come 
first,  then  all  the  machines  in  area  2,  and  so  on. 

These  interconnection  properties  reflect  on  the  characteristics  of 
the  fast  dynamics,  determined  by  the  matrix  of  (3.10).  In  the  remainder  of 
this  section  we  will  study  the  properties  of  B^.  From  equation  (3,10) 

B„ =  A__-L  A, Using  L  we  can  easily  derive  that  the  elements  of  B»  are 

2  22  g  12  g  2 

given  by 


(4.13) 


pq  ak(p) ,q 

nl 

p£  k(p),p’ 

i-j 

st^p 

From  (2.7),  (2.9),  and  (2.10) 


tr  B2  =  tr  A  -  tr  B1 


(4.14) 


it  is  easy  to  see  that  trB^  contains  the  total  sum  of  connections  between 
machines  inside  coherent  areas.  Since  in  the  case  of  slow  coherency 


|tr  B2i  |tr  BJ 


,  IK 

(r  <  y) 


(4.15) 


it  follows  that  the  average  connection  strength  within  areas  is  higher  than 
the  average  interconnection  strength.  To  state  this  simpler:  on  the  average. 
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connections  are  stronger  between  machines  inside  the  same  area  than  between 
machines  from  different  areas.  A  consequence  of  this  property  is  that  the 
matrix  is  also  block-diagonally  dominant,  as  is  A.  Its  diagonal  blocks, 
of  size  equal  to  the  cumber  of  follower  machines  in  the  areas,  are  functions 
of  elements  of  A  corresponding  to  the  connections  of  machines  inside  areas, 
see  (4.15),  and  off-diagonal  blocks  are  functions  of  interconnection  elements 
of  A.  Hence  the  diagonal  blocks  have  larger  norms  than  their  off-diagonal 
blocks.  When  the  interconnection  table  is  formed  for  ,  the  same  way  Table  4.2 
is  formed  for  A,  the  entries  in  the  table  are  almost  identical  to  those  for  A. 


4.5.  Sensitivity  of  Area  Boundaries  to  Parameter  Changes 

In  this  section  we  will  show  that  in  the  case  of  slow  coherency  areas 
are  not  only  weakly  coupled  as  established  in  the  previous  section,  but  also 
area  boundaries  are  relatively  insensitive  to  small  changes  of  network  para¬ 
meters.  The  first  step  in  this  study  is  to  evaluate  the  changes  in  the  steady 
state  matrix  A,  for  small  changes  in  network  parameters.  A  similar  problem 
has  been  treated  in  power  system  analysis  within  so-called  security  analysis. 
Security  analysis  is  concerned  with  the  problem  of  fast  computation  of  the 
equilibrium  state  after  some  given  disturbance.  These  methods,  surveyed  in 
[42],  are  numerical  in  nature,  despite  the  fact  they  use  various  forms  of 
sensitivity  relations  to  achieve  necessary  effectiveness.  There  are  two 
characteristics  which  distinguish  the  approach  in  this  section  from  the 
security  analysis  approach.  First,  our  ultimate  goal  is  to  evaluate  the 

change  of  L  after  a  change  in  network  parameters,  and  second,  we  will  use  the 
§ 

analytical,  rather  than  numerical  approach,  with  the  purpose  of  gaining 


I 


insight  into  the  behavior  of  area  boundaries.  For  that  reason  we  will 
consider  only  the  change  of  admittance  elements  (which  may  model  a  line 
outage,  for  example). 

Outage  of  a  line  between  nodes  i  and  j  affects  matrix  A  in  two 
ways:  first,  it  directly  changes  admittance  y^,y^,y^,y  ,  and  second,  it 

changes  the  equilibrium  angles  6^,i = 1,2, . . . ,n.  Angle  changes  can  be 
analyzed  from  the  DC  load  flow  equations  [52] 


P  =  B  •  9 


(4.16) 


P  -  n-1  vector  of  active  power  injections  into  the  nodes 
of  the  system 

B  -  admittance  matrix, reactive  part 

0  -  n-1  vector  of  angle  differences  6^-6^,  i  =  1,2, . .  .  ,n-l, 
After  a  line  outage  B  changes  to 


b'  =  b  +  ab 


(4.17) 


where 


AB  =  A  b.  .  e  .  e .  . 

ij  ej  ij 


(4.18) 


Ab^  -  the  change  in  admittance  of  the  line  (i,j) 

T 

e^.  -  row  vector  with  all  zero  elements  except  the  i-th,  which 
is  1  and  the  j-th  which  is  -1. 


For  small  Ab„  the  change  in  equilibrium  angles  can  be  found  as  [42] 

A0  =  -B-1AB6. 


(4.19) 


Using  the  symmetry  of  B  a  bound  for  the  change  in  9  can  be  established  as 


A0II  <  Ab,  .(5.-6.)  --  -Ty  ■ . 
ij  i  J  min | A (B) | 


(4.20) 


t 
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This  expression  shows  several  facts:  (i)  in  systems  working  close  to  their 
static  stability  boundaries,  even  small  changes  in  parameters  can  cause  big 
changes  in  equilibrium  angles,  (ii)  parameter  changes  in  higher  loaded  lines, 
and  (iii)  bigger  parameter  changes  can  cause  bigger  equilibrium  changes. 

In  the  rest  of  this  analysis  we  assume  that  parameter  changes  are 
small  so  that  the  linear  analysis  can  still  be  applied  after  the  disturbance. 
Hence  in  determining  sensitivity  of  area  boundaries  we  concentrate  on  the 
effect  of  the  admittance  change  (4.17).  We  want  to  see  the  effect  of  this 
change  on  the  change  in  L  =  L , .  We  approach  the  problem  by  studying  how  the 

8  U  r  i 


slow  subspace  of  A  changes.  From  = L^  it  follows  that  Vo 


I 

=  L  i 

L  gJ 


is  a  basis 


of  the  slow  subspace  of  A.  To  find  a  first  order  variation  of  this  subspace 
corresponding  to  the  disturbed  matrix  A'»M  ^"(B+AB)  we  use  Lemma  2.11  and 
apply  only  one  simultaneous  iteration  to  A'  with  as  the  initial  condition. 

To  make  the  slowest  spectrum  dominant  prior  to  applying  the  simultaneous 
iterations,  we  shift  the  spectrum  of  A,  as  well  as  the  spectrum  of  A',  for 
some  y  to  the  right.  This  is  possible  in  this  case  since  all  eigenvalues  of 
A  are  real  and  in  the  left  hand  side  of  the  complex  plane,  by  assumption.  The 
result  on  the  variation  of  L  is  given  in  the  following  lemma. 

s 

Lemma  4.6:  Suppose  that  for  the  system  (3.11)  cos  terms  in  (3.6)  are 

approximatley  1,  and  that  there  exists  L  =  L,  s.t.  R(L  )  = 0.  Then  for  a  suf- 

g  a  g 

ficiently  small  perturbation  AB  of  the  type  (4.17) 

(i)  if  i  and  j  belong  to  the  same  area,  the  first  order  variation 

of  L,  is  zero 
d 

(ii)  if  i  and  j  belong  to  different  areas,  the  first  order  varia¬ 


tion  of  L^  is  zero  except  for  rows  i  and  j . 


Proof ;  Let  A* =  A+AA  =  M  ^(B+AB).  Due  to  the  linear  convergence  of  the 
simultaneous  iterations  (SI) ,  one  suitably  defined  SI  may  be  used  to  obtain 
a  first  order  variation  of  the  slow  subspace  of  A  when  perturbed  as  A+AA. 

For  that  we  have  to  make  SI  converge  to  the  slow  subspace  of  A*.  This  we  do 
by  using  (2.40)  and  A"  =  A'+yI,  p>max|X(A)|,  instead  of  A.  The  ir  “ial  guess 

slow  subspace  of  A.  After  >  SI, 
(B^+pI)  and  the  fact  that  ,+pI)  is 

nonsingular,  we  can  write 


’  I 

/  =  T  ,  whi 

°  LLgJ 


which  is  a  basis  of  the 


V  =  A"V  . 
1  o 


Using  (M  ^B+pI) 


1  1  1 

Lo  =  L 

L  sJ  g 


(4.21) 


(4.22) 


In  case  (i)  ,  e^VQ  =  0  because  i  and  j  belong  to  the  same  coherent  area.  In 
case  (ii) ,  it  can  be  seen  from  (4.22)  that  the  conclusion  follows  directly. 
Hence,  the  lemma  indicates  that  for  small  parameter  changes  of  the  type  (4.17) 
all  machines  remain  in  their  predisturbance  areas,  except  possibly  the 
terminal  machines  i  and  j .  An  easy  extension  of  the  above  argument  can  be 
used  to  show  that  if  AB  is  modeled  to  account  for  parameter  changes  of  more 
than  one  line  within  one  coherent  area  and  no  changes  in  the  others,  the  only 
nonzero  rows  in  the  matrix  AV  are  those  corresponding  to  the  machines  in  the 
affected  area.  This  means  that  in  the  case  of  a  disturbance  only  machines  in 
the  affected  areas  may  migrate  from  this  area,  and  probably  only  the  boundary 
machines.  Due  to  the  weak  coupling  between  areas  the  situation  as  assumed 


above  may  model  a  line  outage  within  a  coherent  area. 


4 


l 


I 


1 


4 


H 


4 


1 
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A  natural  question  which  arises  after  this  analysis  is  how  can  the 
predisturbance  areas  be  used  in  the  determination  of  areas  after  a  disturbance. 
There  are  several  possibilities. 

(i)  When  the  allowed  computation  time  is  minimal  and  the  distur¬ 
bance  is  small,  the  predisturbance  areas  can  be  used. 

(ii)  The  first  order  variation  formula  (4.22)  can  be  used  to 
decide  how  to  group  borderline  machines.  Any  decision  about  dislocating 
machines  can  be  checked  via  the  criterion  (4.6)  (only  for  the  affected  area 
and  the  machines  in  it) . 

(iii)  Use  the  exact  matrix  A'  (from  a  load  flow  or  static  state 
estimator)  and  Vq  as  initial  guess  to  perform  necessary  number  of  SI.  A 
simplification  in  this  approach  can  be  obtained  if  far  from  disturbance  areas 
are  replaced  by  their  slow  equivalent  derived  in  Section  4.2.  This  is  equi¬ 
valent  to  some  sort  of  condensation  used  in  structural  mechanics  [45]. 

The  following  example  using  the  48  machine  system  confirms  the 

expected  robustness  of  area  boundaries  with  respect  to  the  parameter  and  load 

changes,  and  hence  indicates  that  area  boundaries  will  have  to  be  changed 

only  infrequently.  The  example  corresponds  to  applying  the  alternative  (iii) 

above  for  finding  areas  of  the  network  with  different  configurations.  Three 

network  configurations  are  considered.  For  each  the  exact  L,  is  computed 

d 

using  the  corresponding  linearized  state  matrix.  The  areas  are  obtained  by 
approximating  this  by  an  as  explained  in  Chapter  3. 

Case  one  is  the  nominal  configuration  of  the  system,  which  has  been 
considered  so  far.  For  this  system,  is  given  in  Table  3.1.  The  corresponding 
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areas  are  given  in  Figure  3.7.  Case  two  is  the  same  network,  but  after 

switching  a  line  in  area  1  [45]  following  a  disturbance  in  the  area.  The 

third  case  is  similar  to  case  two,  but  with  a  line  disconnected  in  area  7. 

Inspection  of  the  elements  of  L,  for  all  the  three  cases  shows  that  none  of 

d 

the  machines  changes  its  area  belonging  except  for  the  machines  1  and  2,  which 
in  case  2,  instead  of  being  grouped  with  the  area  1,  become  a  part  of  area  7. 
The  rows  of  corresponding  to  rhree  characteristic  machines,  machines  1,  2, 
and  9,  are  given  in  Table  4.3  for  the  three  network  configurations.  The  dis¬ 
crimination  factors  for  machines  1,  2,  and  9,  also  given  in  the  table,  are 
considerably  smaller  than  the  average  value  of  DF  (for  case  one) ,  which  is  0.85. 
It  is  interesting  that  the  reference  machines  in  all  three  cases  remain  the 


4.6.  Use  of  Slow  Coherency  in  Nonlinear  System  Analysis 

Weak  coupling  between  areas  in  the  case  of  slow  coherency  and  insen¬ 
sitivity  of  area  boundaries  to  small  parameter  changes  indicate  the  possi¬ 
bility  of  using  coherent  areas  of  the  linearized  model  (3.11)  for  the  nonlinear 
system  analysis.  Existing  coherency  methods  [19]  use,  in  fact,  the  linearized 
model  to  identify  coherent  areas,  and  use  them  for  nonlinear  system 
analysis.  In  this  section  we  will  analyze  aggregation  of  the  nonlinear  model 
using  the  block  diagonalizing  transformation  for  the  linearized  one.  This 
analysis  is  along  the  lines  of  the  approach  in  Section  2.5,  except  that  here 

we  take  explicitly  into  account  the  specific  structure  of  L  and  H  .  Later  in 

§  8 

the  section  we  discuss  different  possibilities  for  approximate  system  analysis 


based  on  the  division  of  the  system  into  slow  coherent  areas, 
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The  nonlinear  model  (3.3)  reduced  to  machine  terminals  is 
rewritten  here  for  convenience. 


6  +  M_1D6  =  M-1(P  -P  (6) ) 
m  e 


(A. 23) 


Similarly  to  the  approach  in  Section  2.5,  we  apply  the  transformation 


5  «  6  +  T 


(4.24) 


using  slow  and  fast  deviations  around  the  equilibrium,  y  and  z,  respectively, 
and  the  block-diagonalizing  transformation  of  the  linearized  model  T  from 
(2.52).  In  component  form  (4.24)  becomes 


where 


6  =  6°  +  y  -  Z  z  ,  iSfl 

iej(i)  ai 


6i  =  5i  +  yk-  z 
1  K  i€J(k)  ak 


(4.25) 


k  =  k(i)  -  representative  machine  for  the  area  containing  i 

J(i)  =  I(i)  -  {i} 

M  .  =  Z  M„  . 
al  tei(i)  £ 

After  some  manipulations  with  the  transformed  variables  in  (4.23),  we  get 


y.+D.y.+  Z  D .  z .  =  P  .  -  P  .,  i=  1,2, . 
i  ai  x  . _  T  / . \  1  J  max  eai 
j€J(i)  J  J 


(4.26) 


D.  D  D  D.  M 

h  +  <ML-s£)yi(k)  +  (st-mL)  1  r;i 

3  j  \  J  \  j  iSJ(k)  ak 


P  .  P 


P  .  P  . 


+  -1;  =  +  /_£!__!!£:) 
M .  Zj  V  M .  >L  '  '  M.  M. 

J  J  K  J  K 


j  £  3 


(4.27) 


where 


P  =  IP 

mai  a.ei(i)  m£ 


D  .  =  ID. 
ai  «ei(i)  1 


p  -  IP.. 
eal  Jfil(i)  ei 


(4.28) 

(4.29) 

(4.30) 


Equation  (4.26)  describes  the  area  motion,  while  equation  (4.27)  describes 
the  intermachine  motion.  For  small  variation  around  the  equilibrium  point, 
these  two  equations  are  decoupled.  We  assume  that  even  for  large  excursions 
around  equilibrium,  variables  y  basically  determine  area  motion  and  variables 
z  the  intermachine  motion.  We  now  analyze  in  more  detail  the  area  motion. 
Writing  power  output  of  area  i,  Pgai  as  the  sum  of  intra-  and  interarea  load 
flow,  we  get  (after  a  few  manipulations)  the  accelerating  power  of  the  area  in 
the  form 


F  . (y,z)  HP  .  -  P  , 
ai  mai  eai 


=  -  I  [  I  Z  b .  sin(<5  +* .  (*))]  cosy 

£=1  jei(i)  s€I(£)  2  J  J  X 

IH 

+  Z  Z  Z  b .  sin  5°. 

1=1  j€l(i)  seKi)  JS  JS 


Z  Z  b . „sin(6  +  y  (z) ) 

j€I(i)  «€I(i) 


r 

-  Z  (  Z  Z  b .  cos (6 . 

*-l  j€I(i)  sSI(£)  JS  JS 

it  i 


+  *.s))s1nylt 


E  P»ai(y-z)-Peal(>’’!!) 


(4.31) 


i 


103 


where 


P  .  -  the  first  three  terms  in  (4.31) 
mai 

Pgai  -  the  fourth  term  in  (4.31) 


Z,  ly-VV 


From  the  linear  system  analysis,  the  following  is  true. 

Proposition  4.7:  Assume  L  and  H  satisfy  (2.13)  and  (2.49),  for  system 

8  8 

(3.11)  .  Then: 


(i) 

w°-°>  ■ 0 

(4.32) 

(ii) 

dzFal(0,Z)  -  0 

(4.33) 

(iii) 

yMi<p-0)  ■  °- 

(4.34) 

Proof :  (i)  and  (iii)  follow  directly  from  (4.31).  The  property  (ii)  follows 

from  the  block  diagonalizing  property  of  the  transformation  T  for  the 
linearized  model. 

Based  on  these  properties  we  make 

Assumption  4.8:  The  terms  'i'(z)  have  a  small  effect  on  y  and  can  be  neglected 
in  (4.26)  and  (4.31). 

With  this  assumption  the  model  for  area  dynamics  can  be  realized 
as  a  machine-impedance  system,  as  explained  in  Section  4.2,  but  with  time 
varying  input  P  . (y,0).  The  other  possibility  exploited  in  [43]  is  to 

Kiel  X 

sacrifice  the  symmetry  of  the  network  for  the  advantage  of  having  the 
constant  mechanical  input.  In  this  model,  accelerating  power  is  written  as 


Fai(y>0)  "  -t‘1Y1tsin(’,il  +  ’u)  +  Pmai 


(4.35) 


jei(i)  sei(*)  Js  J! 


b . „  ■  £  E  b,  cos  S . 

11  j€I(i)  361(4)  Is  3s 


(4.38) 


arc  t 


(4.39) 


In  this  model  power  input  is  constant  and  each  line  requires  one  phase 

shifter,  which  eliminates  admittance  symmetry.  Let  us  now  analyze  the 

intermachine  dynamics.  In  equation  (4.27)  the  following  is  assumed. 

°k 

Assumption  4.9:  The  contribution  of  ( -  — ) • ( • )  is  smaller  than  the  contri- 
D  j  It 

bution  of  rp  z  in  equation  (4.27)  so  that  it  can  be  neglected. 

j  J 

Then,  the  equation  of  intermachine  motion  becomes 


where 


(S'A)  +-i  (6,-5,)  =  P„^-P 


J  k  Mj  J  k 


P  .  P  . 
mi  rak 


md  M. 


md  ed 


(4.40) 


(4.41) 


P  .  P  , 
ei  ek 


ed  M. 


(4.42) 


Equations  (4.40)— (4.42)  hold  for  any  partial  coherency.  However,  in  the  case 
of  slow  coherency  furhter  simplification  is  possible  by  using  the  weak 
coupling  between  areas.  In  this  case  it  can  be  assumed  that  the  fast 
variables  are  excited  (different  from  zero)  only  in  the  area  in  which  the 


disturbance  has  occurred.  The  model  for  fast  dynamics  now  becomes 
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P  .  +  ~  Z  b. .sin (6., -6.,  ) 
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where 


E  b.  , 
£Sl(k)  k£ 

(#k 


sin  6  -  P 

£k  o 


(4.43) 


<3^  -  is  the  set  of  all  follower  machines  in  the  affected  area  k 

and  P^  is  artificially  added  to  satisfy  the  equilibrium 

condition  (6..  =6.,  =0). 

Jk  jk 

Equation  (4.43)  can  be  solved  independently  from  the  equation  for 
the  slow  dynamics.  However,  the  solution  angles  are  with  respect  to  the 
reference  machine  angle,  which  becomes  known  only  after  the  slow  subsystem 
(4.26),  (4.31)  is  solved  and  the  transformation  (4.25)  is  used. 


4.7.  48  Machine  Example 

This  section  starts  with  a  brief  discussion  of  the  results  already 
given  for  the  48  machine  system  in  various  sections  of  Chapters  3  and  4.  In 
addition  we  will  also  use  system  responses  to  illustrate  some  of  the  developed 
theoretical  points. 

In  Section  3  the  48  machine  system  has  been  divided  into  nine  areas, 
based  on  nine  slow  modes  and  the  corresponding  approximated  by  an  L^.  Our 
choice  of  slow  coherency,  rather  than  any  other  coherency  was  motivated  by 
the  weak  coupling  property  of  the  areas,  demonstrated  in  Table  4.2  and 
the  robustness  of  area  boundaries,  demonstrated  in  Table  4.3.  These  pro¬ 
perties  of  slow  coherent  areas,  derived  on  a  linearized  model,  will  be  further 
substantiated  in  this  section  by  using  system  responses  of  the  nonlinear  model. 


First  we  check  the  definition  of  slow  coherency  as  given  by 

equation  (3.29).  For  that  purpose  a  fault  is  applied  in  area  one  followed 

by  switching  off  a  line.  The  definition  of  areas  is  kept  the  same  as  in 

the  predisturbance  linearized  system.  Exact  nonlinear  responses  are  obtained 

and  then  the  slow  coherency  definition  checked  by  forming  the  differences  of 

responses  of  machines  inside  the  areas.  As  expected  the  area  with  the  most 

excited  intermachine  motions  is  the  disturbed  area,  i.e.  the  area  one.  Some  z 

variables  for  this  area  are  given  in  Figure  4.1  (see  also  [45]),  and  their 

slow  components  are  given  in  Figure  4.2.  It  can  clearly  be  seen  that  the 

difference  variables  contain  some  slow  dynamics;  basically  all  the  individual 

machine  responses  agree  in  their  slow  motion.  In  our  definition  of  slow 

coherency  we  tolerate  an  arbitrary  magnitude  of  z  variables.  This  is  one  of 

the  key  differences  with  the  other  coherency  definitions  [19],  which  restrict 

the  magnitude  of  z  also.  In  this  particular  example  for  all  the  machines  of 

area  one  to  be  in  the  same  coherent  area  according  to  the  latter  criterion 

the  tolerance  should  be  (see  Figure  4.1)  about  25°,  but  with  this  tolerance 

it  may  as  well  happen  that  many  areas  collapse  into  one  area.  With  our 

approach  of  identifying  areas  based  on  approximation  of  L,  by  an  L  ,  we  avoid 

d  g 

a  need  for  specifying  the  above  tolerance  criterion,  which  may,  as  indicated 
above,  significantly  affect  the  found  areas. 

The  next  figures  also  illustrate  the  robustness  of  area  boundaries 
to  parameter  changes,  which  for  the  linear  case  was  illustrated  by  Table  4.3. 
For  this  purpose  a  fault  was  applied  in  area  7 ,  followed  by  switching  a  line 
in  the  same  area  [45].  We  now  give,  for  the  sake  of  comparison,  the  responses 
of  machines  in  area  7.  Figure  4.3  contains  machine  angle  differences  and 
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Figure  4.4  contains  their  slow  components.  Again  as  in  the  previous  case 
the  definition  of  slow  coherency  is  satisfied  for  area  7  as  well  as  for  all 
other  areas . 

The  third  group  of  figures  illustrates  results  of  various  approxi¬ 
mations  in  the  simulation  of  the  nonlinear  model. 

Approximation  1  (Al) :  In  this  approximation  the  singular  perturbation  approach, 
similar  to  the  one  in  Chapter  2,  was  applied  to  the  model  (4,26),  (4.27). 

Namely,  the  fast  system  (4.27)  is  assumed  to  be  stable  for  each  value  of  y(t) 
and  infinitely  fast,  so  that  the  set  of  equations  containing  differential 
equations  (4.26)  and  algebraic  equations  (4.27)  obtained  by  setting 
z(t)  =  z(t)  =  0  is  solved.  Since  these  equations  describe  the  slow  dynamics, 
a  large  integration  step  size  can  be  used.  However,  the  number  of  equations 
is  still  n.  From  Figure  4.5a  it  can  be  seen  that  the  individual  machine 
responses  are  well  approximated. 

Approximation  2  (A2) :  The  weak  coupling  between  areas,  which  implies  the  weak 
coupling  between  the  fast  subsystems  in  different  areas,  suggests  using 
the  differential  equations  for  the  area  affected  by  the  disturbance,  only, 
and  algebraic  equations  for  all  the  other  areas.  The  number  of  differential 
equations  equals  r+  the  number  of  the  free  machines  in  the  affected  area. 

Again,  the  total  number  of  equations  equals  n.  Inclusion  of  the  fast 
dynamics  in  the  affected  area  improves  the  accuracy  of  the  responses  over  the 
previous  case,  as  can  be  seen  by  comparing  Figure  4.5a  and  Figure  4.5b.  The 
difference  between  the  exact  response  and  A2  is  only  due  to  the  fast  dynamics 
outside  the  affected  area. 

Approximation  3  (A3) :  This  approximation  differs  from  A2  in  that  outside 
the  affected  area  a  constant  z(t)=z(0  )  is  used.  However,  more 
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Figure  4.5.  Exact  individual  machine  responses  in  area  1  vs.  different  approximations.  Al-algebraic 
equation  for  the  difference  variables;  A2-dif f erential  equations  for  the  difference 
variables  in  area  1  only,  algebraic  elsewhere;  A3-differential  equations  for  the 
difference  variables  in  area  1  and  constant  values  for  the  other  areas. 


appropriate  would  be  to  use  z(t)  =  z(“),  but  this  requires  solving  an  additional 
load  flow.  The  results  given  in  Figure  4.5c  indicate  that  this  approximation 
is  somewhat  worse  than  A2. 

As  a  conclusion  to  this  series  of  experiments,  it  follows  that  the 
approximation  A3  seems  to  be  a  compromise  between  accuracy  and  integration  cost. 
Note  however,  that  these  simulations  cannot  indicate  only  the  instability  due 
to  the  fast  dynamics  outside  the  affected  area  because  in  all  the  approximations 
it  was  assumed  that  the  fast  dynamics,  i.e.,  intermachine  motion  is  stable. 

Due  to  the  weak  coupling  between  areas  this  does  not  appear  to  be  a  serious 
limitation.  Therefore  a  justifiable  approach  to  the  system  stability  would 
be  to  use  the  dynamic  equations  for  the  area  dynamics,  and  the  fast  dynamics 
for  the  affected  area  only. 


4.8,  Numerical  Aspects  of  the  Slow  Coherency  Algorithm 

Slow  coherency  proves  to  be  easier  to  identify  numerically  than 
other  types  of  coherency.  In  view  of  the  grouping  algorithm  of  Section  3.5, 
this  is  to  say  that  it  is  somewhat  easier  to  find  a  basis  for  the  slow  eigen- 
space  than  for  any  other  arbitrary  eigenspace  (excluding  the  fast  eigen- 
space  which  is  physically  unrealistic  to  be  used  for  coherency  definition) . 
This  section  is  devoted  to  the  discussion  of  various  methods  for  computing 
a  basis  for  the  slow  subspace,  with  the  emphasis  on  large  and  very  large 
systems.  By  large  systems  we  mean  systems  in  which  the  computation  of  all 
eigenvalues  and  eigenvectors  is  expensive,  but  their  state  matrix  A  can  still 
be  manipulated  in  the  core  of  an  available  computer.  Very  large  systems  will 


be  those  in  which  the  state  matrix  cannot  be  processed  in  core;  computation 
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of  all  eigenvalues  and  eigenvectors  is  almost  impossible.  Independently  of 
the  system  size  there  are  two  factors  that  make  computation  of  V  easier. 

(i)  Matrix  A  is  similar  to  a  symmetric  matrix  (3.14).  This 
means  that  special  programs  can  be  used  which  are  faster  and  more  economical 
in  memory  use  than  the  programs  for  general  matrices. 

(ii)  Arbitrary  basis  of  the  slow  eigenspace  can  be  used,  rather 
than  the  eigenvector  basis.  In  the  case  of  identical  or  very  close  eigen¬ 
values,  eigenvectors  are  very  sensitive  to  parameter  changes  (numerical  error, 
say),  while  the  space  they  span  is  not.  Therefore  any  basis  of  that  subspace 
can  be  used.  Most  of  the  existing  eigenvector  programs  compute  efficiently 
an  arbitrary  basis  corresponding  to  a  multiple  eigenvalue  and  then  spend  most 
of  the  time  trying  to  compute  exact  eigenvectors. 

For  the  large  systems  there  are  already  available  production 
programs  [46]  which  compute  V.  The  conclusion  [46]  is  that  if,  roughly,  r<^r, 
which  can  be  considered  satisfied  in  coherency  identification,  then  the 
special  purpose  programs  for  partial  eigensolutions  are  superior  to  the 
ones  that  compute  all  eigenvalues  and  eigenvectors.  Following  recommendations 
of  [46],  the  bisection  method  seems  to  be  a  desirable  choice.  In  this  method 
the  frequency  range  defining  coherency,  rather  than  the  number  of  slow  modes, 
should  be  specified. 

Alternatively,  a  very  simple  method  of  Chapter  3,  based  on  Riccati 
iterations  (2.41),  can  be  used.  For  this  method  an  initial  guess  of  reference 
machines  needs  to  be  made.  Good  guesses  result  in  faster  convergence  to  the 
solution;  poor  guesses  may  require  many  more  iterations,  as  is  evident  from 
Figure  4.6.  Curve  a  corresponds  to  a  good  choice  of  reference  machines, 
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Figure  4.6.  Behavior  of  Riccati  iterations  for  the  48  machine  system; 

(a)  a  good  guess  of  reference  machines,  (b)  a  poor 
guess  of  reference  machines. 
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in  fact  to  those  that  are  found  by  the  grouping  algorithm.  Curve  b  corre¬ 
sponds  to  selecting  machines  13  and  47  instead  of  17  and  29  in  the  above 
choice.  In  this  case  machines  13  and  47  both  belong  to  the  same  area  (area  7) 
and  area  8  does  not  have  any  reference  machine.  If  this  were  the  ideal 
coherency  case,  the  algorithm  would  not  have  converged  at  all. 

Computation  of  partial  eigensystems  for  very  large  systems  is 

still  an  area  of  active  research  [47,48,32].  Almost  all  the  methods  of  this 

T 

group,  of  which  we  will  mention  only  a  few,  have  inner  product  a  x  as  a 
basic  in  core  operation.  Methods  with  this  property  are  simultaneous  iteration 

methods  [49,32],  modified  simultaneous  iteration  methods  [48],  Lanchos 
method  [49],  and  modified  inverse  method  [47].  All  of  these  methods  have 
basically  the  linear  rate  of  convergence.  Some  acceleration  is  possible  by 
using  Tchebishev's  transformation  [49].  Also,  all  of  these  methods  can  take 
advantage  of  the  sparsity  in  programming  the  inner-product.  However,  in  the 
model  (3.11),  the  state  matrix  A  is  dense.  We  now  show,  on  the  example  of 
simultaneous  iterations,  that  it  is  possible  to  use  the  original  sparse 
Jacobian  (3.4),  instead  of  the  reduced  one.  The  result  is  given  in  the 
following  lemma. 

Lemma  4.10:  Let  J  be  the  Jacobian  of  the  unreduced  load  flow  system  in  which 
the  first  n  rows  correspond  to  the  generator  nodes  and  the  last  m  rows  to  the 
load  nodes.  Let 

,  J  =  JT,  AG  Rnxn  (4.44) 

and  assume  D  ^  exists.  Then 
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(4.45) 


(4.46) 


(4.47) 


generates  for  sufficiently  large  k  a  basis  of  the  r-dimensional  slow  eigen- 
space  of  M  ^ (A-BD  ^C)  =M  ■'"K,  where  K  is  given  by  (3.8),  in  the  form  of  the 
matrix  of  (4.46). 

Proof :  Writing  (4.45)  in  expanded  form,  solving  for  from  the  second 
equation,  and  substituting  in  the  first  equation,  gives  in  view  of  Lemma  2.11 
the  claimed  result. 

Use  of  the  sparse  Jacobian  (4.44)  rather  than  the  reduced  one  as  in 
(3.11),  can  result  in  substantial  savings  in  computation  time,  as  is  already 
shown  in  related  approaches,  see  for  example  [56], 


4.9.  Summary  of  the  Chapter 

In  this  chapter  we  have  shown  that  the  equivalent  for  the  area 
dynamics  in  the  linearized  model  can  be  realized  by  a  symmetric  network 
connecting  r  generators,  each  representing  one  area.  In  the  nonlinear  model 
the  same  structure  can  be  retained  under  very  mild  assumption,  but  with  time 
varying  mechanical  input,  or  alternatively,  asymmetric  network  and  constant 
input.  It  is  concluded  further  that  for  most  networks  and  loading 
conditions,  the  coherency  is  largely  determined  by  network  parameters 
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and  machine  intertias,  rather  than  by  loading  conditions.  The  most  important 
result  is  that  in  the  case  of  slow  coherency  areas  are  weakly  coupled. 
Furthermore,  in  this  very  case  the  area  boundaries  are  relatively  insensitive 
to  small  parameter  changes.  The  weak  coupling  and  robustness  of  area 
boundaries  are  a  basis  for  using  areas  obtained  in  the  analysis  of  the 
linearized  model  for  the  nonlinear  system  analysis.  Several  approximations 
based  on  the  use  of  slow  coherency  and  singular  perturbations  on  the  nonlinear 
model  are  demonstrated  to  be  accurate  and  efficient.  Numerical  aspects  of  the 
slow  coherency  are  discussed  in  light  of  the  existing  methods  and  research 
trends  in  the  field.  A  sparse  formulation  of  the  numerical  algorithm  for 
computing  slow  basis  is  given. 
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5.  SUGGESTIONS  AND  CONCLUSIONS 

5.1.  Suggestions  for  Further  Research 

In  chapter  4  we  gave  several  uses  of  the  coherency.  First 
we  have  considered  aggregation  of  linear  and  nonlinear  power  system 
models.  Then  by  using  slow  coherency  we  were  able  to  transform  the  mixed 
state  electromechanical  power  system  model  into  the  singularly  perturbed 
form.  Finally,  we  illustrated  several  uses  of  the  model  in  the  singularly 
perturbed  form  as  well  as  the  weak  coupling  between  coherent  areas  in 
the  simulation  of  the  nonlinear  model.  A  simulation  of  the  system  is 
usually  undertaken  for  the  purpose  of  checking  its  stability.  We  now 
want  to  discuss  possible  uses  of  coherency  and  in  particular  slow  coherency 
in  the  direct  stability  analysis  of  power  systems.  In  order  to  do  so 
we  will  first  sketch  the  background  of  the  direct  stability  analysis. 

One  of  the  main  objectives  in  the  stability  studies 

of  power  systems  is  to  determine  the  so  called  critical  clearing  time  tc> 
This  is  the  time  defined  for  each  fault  and  represents  the  maximal 
duration  of  the  fault  for  which  post  disturbance  system  (which  may  have 
different  structure  than  the  original  one)  remains  stable.  Large  t£ 
indicates  that  the  system  is  stronger  (more  robust)  with  respect  to  the 
given  disturbance.  Critical  clearing  time  is  then  used  in  the  design 
of  protective  equipment  (local  controls).  The  problem  of  determining  the 
critical  clearing  time  is,  however,  almost  identical  to  the  problem  of 
determining  the  domain  of  attraction  around  postdisturbance  equilibrium  - 
or  stability  domain  in  short. 

Stability  analysis  of  power  systems  at  the  present  time  is 
basically  conducted  using  simulation  of  system  responses  for  given 
disturbances  [59].  Such  methods  enable  use  of  detailed  models  of  system 


components  (machines  and  transmission  elements)  and  offer  flexibility  in 
using  different  configurations  of  the  same  network.  However  there  are 
some  disadvantages:  (a)  computation  time  is  large,  becoming  prohibitive 
for  very  large  systems,  (b)  experience  needed  in  interpreting  responses 
in  order  to  conclude  instability  and  (c)  no  advantage  from  previous 
simulations  can  be  gained  when  analyzing  the  same  network  for  different 
fault  locations .  An  alternative  is  to  perform  a  direct  stability  analysis 
based  on  Liapunov  theory  [53].  The  electromechanical  model  (3.2),  (3.3), 
in  which  the  load  nodes  are  eliminated  and  the  admittances  between  the 
remaining  generator  busses  are  suitably  altered  as  discussed  in 
section  3.2,  is  typically  used  in  such  studies.  By  further  manipulations 
[53]  the  model  can  be  written  in  the  form 

x  =  Ax  +  B  f(y) 

y  =  Cx  (5.1) 


where  x  =  (u^,^,.. 


.  >'V6rV62'6N’  *  ’  ’  ,6N-1  N  ) : 


f .  (y)  =  V  V  B  [sin(y.  +6°  )  -  sin  6°  ] 

i  p  q  pq  i  pq  pq 


(5.2) 


The  relation  between  indexes  i,  p  and  q  and  the  form  of  matrices  A,  B 
and  C  are  given  in  [53].  For  the  model  (5.1)  of  the  postdisturbance 
network  the  form  of  a  Liapunov  function 

y 

V (x)  =  XTPx  +  2 j*  fT(u)Qdu  =  vc  +  vp  (5.3) 


and  the  necessary  definitnes  conditions  are  given  in  [62].  From  (5.3) 

jM  0  ! 

the  classical  transient  energy  function  [61]  is  obtained  with  P  =  q  q 


and  Q  =  I.  Having  selected  a  Liapunov  function,  the  next  step  is  to 
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find  the  critical  value  V(x)  =  such  that  if  V(x)  <  for  some  x, 
then  the  system  is  said  to  be  stable.  In  [60]  it  is  shown  that  equals 
the  minimum  potential  energy  over  all  the  unstable  equilibrium  points 
of  the  system  (5.1),  i.e. 

V  =  inf  V (x)  (5.4) 

°  x€U 

where  U  is  the  set  of  all  the  unstable  equilibrium  points  for  the  system 
(5.1).  With  known  Vc  computation  of  the  critical  clearing  time  for 
each  disturbance  proceeds  as  follows.  When  the  system  is  subjected  to 
the  disturbance,  its  response  is  simulated  and  the  function  V(x)  of 
(5.3)  is  evaluated  along  the  trajectory.  The  time  at  which  V(x)  =  is 
the  critical  clearing  time.  From  this  description  it  is  clear  that  the 
direct  methods  are  in  fact  the  combination  of  a  simulation  and  a  direct 
stability  analysis.  An  inherent  difficulty  in  the  direct  methods,  however, 
is  computation  of  because  the  number  of  unstable  equilibrium  points 
of  (5.1)  is  0(2n),  where  n  is  the  number  of  generators  in  the  network, 
and  to  find  each  of  this  points  requires  solving  the  set  of  (n-1) 
nonlinear  (load  flow)  equations.  Therefore  the  utility  of  the  method 
critically  depends  on  the  way  of  computing  Vc-,  its  effectiveness  is  to 
be  measured  against  possible  savings  in  computer  time  over  the  direct 
simulation.  There  are  also  other  problems  summerized  below  which  have 
resulted  in  certain  scepticism  from  power  industry  toward  the  application 
of  these  methods  [Disc,  on  61]. 

(a)  The  methods  have  proven  useful  for  small  order  systems  (several 
machines,  say);  for  one  machine  infinite  bus  system  exact  stability 
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boundary  can  be  obtained. 


(b)  Depending  on  the  fault  location  they  can  give  overly  conservative 
results  when  applied  to  large  systems  [56];  if  parts  of  the  same  system 
are  aggregated.,  better  results  can  be  obtained  [56]. 

(c)  Most  of  the  results  for  large  systems  are  reported  for  the  electro¬ 
mechanical  model  and  the  line  conductances  neglected.  Inclusion  of  the 
more  detailed  mechine  models  and  the  conductances  is  followed  by  the 
difficulty  in  finding  an  "appropriate"  Liapunov  function. 

Despite  the  shortcomings,  the  need  for  direct  stability  methods  which 

would  complement  simulation,  exists.  The  utility  of  the  approach  is 

seen  in  the  planning  stage  or  on-line  security  assessment  where  it  would 

be  used  to  select  some  critical  disturbances  for  further  study.  Here 

the  initial  effort  spent  on  computing  would  be  well  compensated  by  the 

repeated  use  of  the  same  V  for  different  faults  resulting  in  the  same 

c 

postdisturbance  network.  These  methods  also  offer  a  basis  for  studying 

the  effect  of  parameter  changes  on  systems  stability.  In  response  to 

this  need,  and  in  an  attempt  to  overcome  (some  of)  the  shortcomings, 

there  have  recently  appeared  several  major  contributions  to  the  problem 

of  direct  transient  stability  studies. 

First,  in  [63]  it  was  shown  that  the  potential  function 

(the  integral  part  of  V(x)  in  (5.3))  is  smooth  around  equilibrium  points. 

Using  this  property  of  the  potential  function  and  generalizing  the  case 

o 

of  c-ne  machine  infinite  bus  system,  in  which  6  is  the  stable  equilibrium, 
then  the  smaller  of  (+n-26°)is  the  unstable  equilibrium  with  the 
smallest  potential  energy.  Approximations  of  the  unstable  equilibrium 
points  of  the  following  type  have  been  suggested  by  several  authors 
(See  [53]  for  an  overview) 
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(5.5) 


where  superscript  "o"  denotes  equilibrium  values,  and  uk  is  an  approxi¬ 
mation  of  the  angle  corresponding  to  the  k-th  machine  becoming  unstable. 
Variants  exist  [53]  to  account  for  several  machines  becoming  unstable. 

Now  the  search  for  in  (5.4)  can  be  carried  out  over  the  set  u  whose 
elements  are  of  the  form  (5.5).  This  set  may  still  be  large  (it  has  n 
elements  if  only  one  machine  at  a  time  is  considered  to  become  unstable, 
but  combinations  often  have  to  be  considered  as  well)  so  that  selecting 
the  right  element  remains  a  serious  task.  It  is  interesting  to  note, 
however,  that  when  the  separation  of  a  system  into  groups  under  a  given 
distrubance  is  known,  then  the  Vc  computed  as  the  potential  energy  at  the 
points  of  the  type  (5.5)  corresponding  to  these  groups  gives  favorable 
results.  A  major  step  forward  in  finding  for  a  given  disturbance  has 
been  made  in  [54,55]. 

The  main  idea  is  that  Vc  is  the  potential  energy  at  the  unstable 
equilibrium  point  which  is  in  the  direction  of  the  systems  motion  (the 
motion  is  considered  in  the  space  of  angles  6).  The  direction  is  determined 
at  the  time  when  the  kinetic  energy  is  maximal.  This  implies  that  the 
system  continues  to  move  basically  in  the  direction  established  during 
the  fault,  even  after  the  fault  is  cleared.  As  a  consequence  the  method 
proves  to  be  very  effective  for  the  so  called  first  swing  instability, 
which  is  when  6  grows  without  changing  signs.  On  the  other  hand  it  is 
often  inaccurate  for  multiswing  instability,  in  which  case  the  system 


oscillates  several  times  before  becoming  unstable.  In  this  case  the 
assumption  of  unidirectional  motion  is  not  satisfied,  and  a  different 
unstable  equilibrium  has  to  be  found  in  order  to  compute  appropriate  V  . 

A  group  of  methods  with  sound  analytical  background,  but  less  satisfying 
practical  results  so  far,  is  centered  around  the  use  of  vector  Lyapunov 
function  approach  [53,57,43].  The  essence  of  the  method  is  to  [58]: 

(a)  decompose  the  system  into  suitable  subsystems,  (b)  find  stability 
domain  of  subsystems  ignoring  the  interconnections,  (c)  impose  (linear) 
bound  on  interconnections  and  (d)  find  the  stability  region  of  the 
overall  system  [57]. 

For  the  methods  to  be  successfully  applied  it  is  essential 

that  the  subsystems  are  weakly  coupled  [57,58].  The  main  reasons  for  the 

conservativeness  are  :  (a)  in  bounding  interconnections  it  is  assumed 

that  they  all  act  simultaneously  in  the  way  which  is  the  worst  for  the 

systems  stability,  and  (b)  nonlinear  interconnection  terms  of  the  form 

s  =  [sin(x-ta)  -  sina  ]  are  usually  bounded  by  ex  5  S  ^  cos  a  *  x, 

where  £  has  to  satisfy  two  contradictory  requirements:  for  a  large 

region  of  stability  to  be  obtained  it  has  to  be  as  small  as  possible, 

*  T 

but  for  the  matrix  A(£),  from  V  =  A(£)V  where  V  is  the  vector  of 
subsystem  Liapunov  functions,  to  have  all  the  eigenvalue  positive,  £  has 
to  be  as  large  as  possible.  It  may  happen,  therefore,  that  for  a  given 
£,  which  is  a  design  parameter,  suitable  A(£)  does  not  exist;  then  £ 
has  to  be  increased,  resulting  ultimately  in  the  reduction  of  the  region 
of  stability.  Ihe  methods,  however,  can  be  very  attractive  when  the 
decomposition  of  the  system  into  weakly  coupled  subsystems  can  be  made. 


With  this  somewhat  longer  introduction  into  the  complex  and  unresolved 
problem,  we  are  now  ready  to  consider  coherency  in  the  framework  of  the 
direct  stability  analysis.  Our  consideration  is  basically  oriented  toward 
listing  problems  for  further  research. 

For  the  methods  based  on  vector  Liapunov  function  approach  the 
decomposition  of  the  system  into  areas  based  on  slow  coherency  seems 
to  be  the  best  choice  because  of  the  weak  coupling  between  areas  established 
in  chapter  4.  If  the  problem  consists  of  finding  the  "slow"  stability 
i.e.  the  stability  of  the  interarea  motion,  assuming  that  the  fast 
intermachine  oscillations  within  areas  are  stable,  then  the  vector 
Lyapunov  function  approaches  appear  to  be  applicable.  However,  when  the 
stability  of  the  overall  system,  including  the  intermachine  oscillations, 
is  of  interest  (which  is  the  most  usual  case),  then  there  are  several 
options.  One  is  to  study  separately  the  fast  and  the  slow  stability  using 
models  (4.26)  and  4.43).  For  analyzing  the  fast  stability,  a  method 
along  the  lines  of  [54,55]  will  have  to  be  used.  The  other  way  is  to 
aggregate  all  the  areas  outside  the  area  in  which  a  distrubance  is  applied 
and  study  the  whole  composite  system  [43] . 

In  summary,  the  coherency  can  help  in  three  ways.  First, 
coherency  based  aggregation  reduces  the  order  of  the  system.  Second, 
when  slow  coherency  is  used,  the  areas  are  weakly  coupled,  and  third,  the 
machines  in  the  coherent  areas  have  a  sort  of  symmetry  (homogenity)  in 
coupling  (4.6).  The  last  property  may  imply  that  the  criterion  (5.4), 
with  U  containing  unstable  points  (5.5)  corresponding  to  the  generators 
in  the  area  affected  by  a  disturbance,  may  not  be  too  conservative. 


Among  other  applications  of  the  coherency  decomposition  and 
possible  extensions  of  the  work  we  mention  the  following. 

The  electromechanical  model  without  damping  was  used  in  the 
previous  analysis  under  the  assumption  that  small  damping  does  not 
affect  frequency  of  the  system.  However  the  effects  of  the  more  complete 
machine  models  on  the  identity  of  coherent  areas  should  be  investigated. 

Numerical  aspects  of  the  simulation  algorithm  for  weakly 
nonlinear  systems,  the  algorithm  for  identifying  coherent  areas  and  the 
use  of  the  decomposition  for  analysis  and  design  still  deserve  a  lot  of 
attention. 

When  the  slow  coherency  is  used  for  system  decomposition,  the 
areas  are  weakly  coupled.  It  is  felt  that  this  property  and  a  similar 
approach  can  be  used  for  the  decomposition  of  a  power  system  for  use  in 
steady  state  (load  flow)  type  studies. 

5.2.  Conclusions 

In  this  thesis  the  properties  and  some  uses  of  coherency  in 
electromechanical  models  of  power  systems  have  been  considered. 

It  has  been  shown  that  the  coherency  can  be  treated  in  the  framework  of 
the  spectrum  separation  problem.  This  property  was  essential  in  defining 
the  new  efficient  algorithm  for  the  decomposition  of  power  systems  into 
coherent  areas. 

When  the  coherency  is  defined  as  the  so  called  slow  coherency, 
then  the  areas  are  weakly  coupled  and  insensitive  to  small  parameter 
changes.  This  is  a  justification  for  using  linearized  model  for  defining 
areas,  and  then  using  those  areas  in  the  nonlinear  model. 
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For  the  class  of  r-decoraposable  power  systems  with  respect  to 
the  slow  modes  the  subsystem  whose  states  are  so  called  center  of 
intertia  variables  is  the  true  slow  subsystem  and  th«  subsystem  with  inter¬ 
machine  variables  within  areas  as  states  is  the  true  fast  subsystem. 
Applying  the  area  and  intermachine  variable  transformation  to  near 
r-decomposable  systems  transforms  the  mixed  state  electro¬ 
mechanical  model  into  the  singularly  perturbed  form.  Use  of  the  standard 
singular  perturbation  approaches  to  the  simulation  of  such  models  was 
illustrated  on  several  examples.  Analysis  of  the  reduced  order  models 
has  revealed  the  sources  of  approximations  used  in  going  from  the  exact 
subsystems  to  the  commonly  used  aggregated  models,  both  in  the  linear 
and  the  nonlinear  case. 

The  numerical  algorithm  for  identifying  coherent  areas  is 
designed  with  the  idea  to  be  applicable  for  large  power  systems,  and  hence 
it  consists  of  well  conditioned  steps  and  modest  use  of  computer  time 
and  memory.  Together  with  the  weak  coupling  property  of  the  resulting 
areas  the  met  jd  offers  a  useful  way  of  decomposing  a  power  system  into 
subsystems  which  can  then  be  used  in  various  applications. 
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APPENDIX  A.  ONE  MACHINE  INFINITE  BUS  MODEL  [34] 

Notation 

AVR  -  automatic  voltage  regulator 

E , ,V  ,R.  -  voltage  regulator  states 
fd  r  f 

t  f 

e  ,e,  -  flux  decay  states 
q  d 

£1,  6  -  swing  model  states 

Subscripts 

d  -  direct  axis  quantities 
q  -  quadrature  axis  quantities 
f  -  main  field  quantities 


eq.Gd 


V, 


1 — U 


Infinite  Bus 


fP-6192 


Figure  A.l.  One  machine  infinity  bus  system. 
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hr. 
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Full  Model 


'  -  [1  +  (x.-x’)Y] 


D 


q  t 


+  ^-Efd> 


-  (xd~x')YVisin6  +  EfdJ 


=  jr~  l(x  -x')  YV£  cos6  -  [1  +  (x  -x')Y]  e^} 

qo 

5  =  377  (1  -  1) 


1  ,  Pin 


^  =  2H  ^  ~q~  -YV^CeqCosS  +  edsin6)] 


£  -i  {  - 


fd  -  IE  '  -  +  W’  Etd  +  V 


kf 


V8-Y[KA(Rf  "V"  "V=+V^f)  -  V 


vc  -  [(l-n'Y)2  (eq2  +  ej2)  +  2 (l-x'Y) (x'YV. ) <ejcos6 )  +  (x'YV  ) 


SE(Efd> 


A  exp  [ (B  )(Er,)] 
sat  ^  1 v  sat''  fd'J 


Y 


1 

x'+xg 


2]  1/2 


134 


For  the  example  of  chapter  2  the  parameters  have  the  following  values 


H 

= 

5.0  sec 

ta 

= 

0.06  sec 

D 

- 

2.0  pu 

te 

= 

0.5  sec 

Xd 

= 

1.2  pu 

tf 

= 

1.0  sec 

X 

q 

- 

1.0  pu 

ka 

= 

25 

x' 

= 

0.25  pu 

“e 

= 

-0.0445 

X 

e 

S 

0.25  pu 

Kf 

* 

0.16 

T' 

— 

5.0  sec 

A 

— 

0.001123 

do 

sat 

T' 

= 

0.50  sec 

B 

as 

0.3043 

sat 

The  operating  point  is  defined  by 


Vi  = 

1.  +  jo. 

p.u. 

Vt  =* 

1*  +  jo. 2 

p.u. 

Sg  = 

0.8  +  jo. 

1608  p.u. 

VRef 

1.  p.u. 

With  these  parameters  linearization  of 
state  matrix 

the  model 

gives  the 

following 

-0.58 

0 

0 

-0.269 

0 

0.2 

0 

0 

-1.0 

0 

0 

0 

0.16 

0 

0 

0 

-5.0 

2.12 

0 

0 

0 

0 

0 

0 

0 

377 

0 

0 

-0.141 

0 

0.141 

-0.2 

-0.28 

0 

0 

0 

0 

0 

0 

0 

0.0838 

2.0 

173 

417 

-116 

40.9 

0 

-66.7 

-16.7 

1 


E 
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APPENDIX  B.  MATRIX  A  OF  6 = A6  FOR  THE  48  MACHINE  SYSTEM 


D 


1 

2 

1 

3 

.36 

4 

5 

6 

7 

. »■ . v . "* — 

8 

1 . 

-0.  150 

0. 037 

0.  010 

0. 004 

0.  003 

0.  012 

0. 009 

2. 

0. 032 

-0.  148 

0.  013 

0. 005 

0. 003 

0.  014 

0. 01 1 

0.  010 

3. 

0. 01 1 

0. 016 

-0. 202 

0.  039 

0. 025 

0.033 

0. 026 

0. 014 

4. 

0.  020 

0. 028 

0.  166 

-0. 583 

0.  172 

0.057 

0. 045 

0.025 

5. 

0. 004 

0. 005 

0.035 

0.  061 

-0. 142 

0. 01  1 

0. 005 

6. 

0. 010 

0. 014 

0. 026 

0.  010 

0.  007 

-0.  172 

0.013 

7. 

0. 010 

0. 014 

0. 027 

0.  010 

0. 007 

0.  080 

-0.196 

0.013 

8. 

0. 006 

0.003 

0. 009 

0.  003 

0.003 

0. 01 1 

0. 009 

-0.089 

9. 

0. 014 

0.  017 

0. 016 

0. 006 

0. 004 

0. 018 

0. 014 

0.034 

10. 

0.007 

0.  007 

0.004 

0.  002 

0. 001 

0.  005 

0. 004 

0.006  ' 

1 1 . 

0. 001 

0. 001 

0. 001 

0.  000 

0.  000 

0. 001 

0. 001 

0.001  • 

12. 

0.  002 

0.  002 

0.001 

0. 000 

0.000 

0. 001 

0. 001 

0. 002 

13. 

0. 002 

0. 002 

0. 001 

0.  000 

0. 000 

0.  002 

0. 001 

0.002 

14. 

0. 003 

0. 004 

0.  002 

0.  001 

0.  001 

0.  002 

0. 002 

0.003  , 

15. 

0.001 

0.  001 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0.001 

16. 

0. 000 

0. 000 

0.000 

0.  000 

0.  000 

0. 000 

0. 000 

0.000 

17. 

-0. 000 

-0. 000 

-0. 000 

-0.000 

-0. 000 

-0. 000 

-0. 000 

0.000  • 

18. 

0.  000 

0.  000 

0. 000 

0.  000 

0.  000 

0.  000 

0. 000 

0.000 

19. 

0. 000 

0.  001 

0.  000 

0.  000 

0. 000 

0. 000 

0. 000 

0. 000 

20. 

0.  000 

0. 000 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0.000 

21  . 

0. 000 

0.  000 

0. 000 

0.  000 

0.  000 

0.000 

0. 000 

0.000  . 

22. 

0.  000 

0. 000 

0.  000 

0.  000 

0. 000 

0. 000 

0. 000 

0.000 

23. 

0.  001 

0.001 

0.  001 

0.  000 

0.  000 

0.  001 

0.  001 

0 . 00 1  ; 

24. 

0. 003 

0. 004 

0.  002 

0. 001 

0.  001 

0. 002 

0.  002 

0.003 

25. 

0. 003 

0. 003 

0. 002 

0. 001 

0. 001 

0. 002 

0. 002 

0. 003 

26. 

0.  002 

0.  002 

0.001 

0.  000 

0.  000 

0.001 

0. 001 

0.002  , 

27. 

0.  005 

0. 005 

0. 002 

0.  001 

0. 001 

0.  003 

0.  002 

0.002  1 

28. 

0.  000 

0.  000 

0. 000 

0.  000 

0.  000 

0.  000 

0.  000 

0.000 

29. 

0.  000 

0. 000 

0. 000 

0. 000 

0.000 

0. 000 

0. 000 

0.000 

30. 

0  000 

0.  000 

0.  000 

0.  000 

0.000 

0.  000 

0.  000 

0.000 

31  . 

0. 000 

0. 000 

0. 000 

0.  000 

0.  000 

0. 000 

0.  000 

0.000  ■ 

32. 

0. 000 

0. 000 

0. 000 

0. 000 

0. 000 

0. 000 

0. 000 

0.  000 

33- 

0. 000 

0. 000 

0. 000 

0. 000 

0.  000 

0.  000 

0. 000 

0.000  - 

34. 

0.  000 

0.000 

0. 000 

0.  000 

0.  000 

0. 000 

0.  000 

0.000  ■ 

35. 

0. 000 

0. 000 

0. 000 

0.  000 

0.  000 

0. 000 

0. 000 

0. 000 

36. 

0. 000 

0.  000 

0. 000 

0. 000 

0. 000 

0.000 

0.  000 

0.000 

37. 

0. 000 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0.000 

0.000  j 

33. 

0. 000 

0. 000 

0. 000 

0.  000 

0. 000 

0.  000 

0.  000 

0.000  1 

39. 

0. 000 

0. 000 

0. 000 

0.  000 

0.  000 

0.  000 

0. 000 

0.  000 

40. 

0.  000 

0.  000 

0. 000 

0. 000 

0. 000 

0. 000 

0. 000 

0. 000 

41  . 

0. 000 

0.000 

0.000 

0.  000 

0. 000 

0.  000 

0. 000 

0.000 

42. 

0. 000 

0.000 

0.000 

-0.000 

0. 000 

0. 000 

0. 000 

0. 000 

43. 

0. 000 

0. 000 

0. 000 

0.  000 

0.  000 

0.  000 

0. 000 

0.  000 

44. 

0. 000 

0. 000 

0.  000 

0.  000 

0. 000 

0. 000 

0.  000 

0.000 

45. 

0. 000 

0.  001 

0.  000 

0.  000 

0. 000 

0. 000 

0.  000 

0. 000 

46. 

0. 000 

0.000 

0. 000 

0.  000 

0. 000 

o.  ooo 

0.  000 

0.000  ■ 

47. 

0. 001 

0.  001 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0.001 

43. 

0.  000 

0.  000 

0.  000 

0. 000 

0. 000 

0.  000 

0.  000 

0.  000 

r"  - 

kr. 

137 
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9 

10 

11 

12 

13 

14 

15 

16 

1  . 

0. 01 1 

0.  002 

0. 001 

0.001 

0.  003 

0.  001 

0.  001 

0. 001 

2. 

0. 012 

0.002 

0.  001 

0.  001 

0. 002 

0.  001 

0. 001 

0.  001 

D 

3. 

0. 013 

0. 001 

0.  000 

0.001 

0.  002 
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0.  001 

0.  000 

4. 

0.023 

0.003 

0.  001 
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0. 003 

0. 002 

0.  001 

0.  001 

5. 

0. 004 

0. 000 

0.  000 

0. 000 

0.  001 

0.  000 

0. 000 

0. 000 

6. 

0. 012 

0. 001 

0.  000 

0.001 

0.  002 

0.  001 

0.  001 

0. 000 

7. 

0.012 

0. 001 

0. 000 

0.001 

0. 002 

0. 001 

0. 001 

0. 000 

8. 

0. 020 

0.  001 

0.  000 

0.001 

0. 002 

0.  001 

0.  001 

0.  000 

m 

9. 

-0. 198 

0. 006 

0.002 

0. 002 

0.  006 

0. 004 

0.  002 

0.  001 

10. 

0. 01 1 

-0.259 

0.  009 

0. 012 

0. 024 

0.  020 

0. 009 

0. 005 

1 1 . 

0. 002 

0. 007 

-0.185 

0.  084 

0. 010 

0. 006 

0. 003 

0.  002 

12. 

0. 003 

0. 010 

0. 083 

-0. 228 

0.  013 

0.  009 

0.  005 

0.  003 

13. 

0. 004 

0. 009 

0. 004 

0.  006 

-0.  172 

0.  027 

0. 015 

0. 009 

14.- 

0. 006 

0. 017 

0. 007 

0. 009 

0.055 

-0. 289 

0. 025 

0.  012 

15. 

0. 001 

0. 002 

0. 001 

0.  002 

0. 012 

0.  009 

-0.  144 

0. 01 1 

16. 

0. 001 

0.  001 

0. 001 

0.  001 

0.  006 

0.  004 

0.  009 

-0.158 

17. 

0. 000 

-0. 002 

0. 000 

0.000 

0. 004 

-0. 003 

0.  007 

0.  020 

18. 

0. 001 

0.  001 

0.  001 

0.  001 

0.  006 

0.  003 

0.  014 

0 . 02  3 

19. 

0. 001 

0. 002 

0. 001 

0. 001 

0. 009 

0. 006 

0. 015 

0 . 04 1 

20. 

0.  001 

0.001 

0. 001 

0.  001 

0. 007 

0. 004 

0.  015 

0. 029 

21  . 

0. 000 

0.  001 

0. 000 

0.  001 

0.  004 

0. 002 

0.  007 

0. 016 

22. 

0.  001 

0.  001 

0. 001 

0.  001 

0.  006 

0. 004 

0. 01  1 

0.024 

23. 

0. 002 

0.  004 

0. 002 

0.  002 

0. 014 

0.013 

0.  069 

0. 015 

24. 

0. 006 

0.  010 

o.  003 

0.  004 

0.  013 

0.  01  3 

0.  018 

0. 01  1 

25. 

0.  005 

0.  009 

0.  003 

0.  004 

0. 01  1 

0. 01  1 

0.  008 

0. 005 

E 

26. 

0. 003 

0. 008 

0. 003 

0.  003 

0. 015 

0.032 

0.  009 

0. 005 

27. 

0.  004 

0. 005 

0.  001 

0.  001 

0.  003 

0.  003 

0. 001 

0. 001 

28. 

0. 000 

0. 000 

0.  000 

0.  000 

0. 000 

0.  000 

0.000 

0. 000 

29. 

0. 000 

0. 000 

0.  000 

0. 000 

0. 000 

0. 000 

0. 000 

0. 000 

30. 

0. 000 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0. 000 

0.  000 

31  . 

0.  coo 

0. 001 

0.001 

0. 001 

0. 005 

0.  003 

0.  003 

0.02  4 

32. 

0. 000 

0. 000 

0.  000 

0. 000 

0. 000 

0. 000 

0. 001 

0.002 

33. 

0. 00 0 

0,  000 

0.  001 

0. 001 

0.  001 

0.  001 

0.  002 

0.  005 

34. 

0. 000 

0.000 

0. 000 

0. 000 

0.  000 

0.  000 

0. 000 

0. 001 

35. 

0. 000 

0. 000 

0.  000 

0. 000 

0.  000 

0.  000 

0. 000 

0 . 00 1 

36. 

0.  000 

0.  000 

0.  005 

0.  008 

0.  001 

0.  000 

0.  001 

0.  001 

37. 

0.  000 

0. 000 

0. 000 

0. 000 

0.  000 

0. 000 

0. 000 

0.  001 

38. 

0. 000 

0. 000 

0.  000 

0. 000 

0.  000 

0. 000 

0.  000 

0. 000 

39. 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

0.  000 

40. 

0. 000 

0. 000 

0. 000 

0. 000 

0. 000 

0.  000 

0. 000 

0. 000 

41 . 

0. 000 

0. 000 

0. 000 

0.000 

0.  000 

0.  000 

0. 000 

0. 000 

42. 

0.  000 

0. 000 

0.  000 

C.  000 
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0. 000 

0. 000 

41 . 
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0. 000 

0. 000 

0.  000 

« 
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0. 000 
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0. 000 
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0. 002 
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0.  002 

0. 00  2 

46. 

0. 000 

0. 000 

0. 000 

0. 000 

0. 00  0 

0.  000 

0.  000 

0. 000 

47. 

0. 001 

0. 002 

0.  001 

0.  001 

0.006 

0. 004 

0. 010 

3.  020 

t 

L 

43. 

0. 000 

0. 000 

0. 000 

0.  000 

0. 000 

0.  000 

0.  00  0 

0.  00  ^ 

. - - . — J 
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0.  000 
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0.  000 

0. 000 
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0.  000 

0. 000 

0. 000 

0. 000 

0. 000 

0.  000 

0.  000 

0.  000 
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0.  000 

0. 000 

0. 000 

0.  001 

0. 000 
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0. 002 

0. 001 
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0.  001 
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0. 001 

0. 007 

0.  004 

0. 002 
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0. 003 
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0.  001 
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0.  002 
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0.  001 

0.  001 

0. 001 

0.  004 

0.  002 
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0.  001 

0.  001 

0. 003 

0. 012 

0.006 

0. 004 

ShwhCh 

lUi  ieB 

0. 016 

0. 00 

8 

0. 005 

0.004 

0. 002 

0. 002 

0. 023 

0.  00 

9 

0. 007 

0.  004 

0. 002 

0. 017 

0. 001 

filfi  ifH 

0. 030 

0. 01 

7 

0.  010 

0. 007 

0. 004 

0.002 
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0. 000 
0.  000 
0. 000 
0.  000 
0. 009 
0. 001 
0. 002 
0. 000 
0. 001 
0. 001 
0. 000 
0. 000 
0. 000 
0. 000 
0. 000 
0. 000 
0. 000 
0. 000 
0. 001 
0.000 

0. 000 


0. 006 
0.  00  1 
0.  000 
0.  000 
0. 000 
0. 028 
0.  002 
0. 006 
0. 001 
0. 001 
0. 001 
0. 001 
0.  001 
0. 000 
0. 000 
0. 000 
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